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initial mesh tandem blade gap region 

7.94 Allison tandem blade cascade 70 % span from shroud 1 st adapted mesh: 

9,367 nodes and 18,186 elements 
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ABSTRACT 


Jiang, Yi-Tsann. Ph.D., Purdue University, May 1993. Development of an Un- 
structured Solution Adaptive Method for the Quasi-Three- Dimensional Euler and 
Navier-Stokes Equations. Major Professor: Dr. William J. Usab, Jr. 

A general solution adaptive scheme based on a remeshing technique is devel- 
oped for solving the two-dimensional and quasi-three-dimensional Euler and Favre- 
averaged Navier-Stokes equations. The numerical scheme is formulated on an un- 
structured triangular mesh utilizing an edge-based pointer system which defines the 
edge connectivity of the mesh structure. Jameson’s four-stage hybrid Runge-Kutta 
scheme is used to march the solution in time. The convergence rate is enhanced 
through the use of local time stepping and implicit residual averaging. As the 
solution evolves, the mesh is regenerated adaptively using flow field information. 
Mesh adaptation parameters are evaluated such that an estimated local numer- 
ical error is equally distributed over the whole domain. For inviscid flows, the 

4 

present approach generates a complete unstructured triangular mesh using the ad- 
vancing front method. For turbulent flows, the approach combines a local highly 
stretched structured triangular mesh in the boundary layer region with an un- 
structured mesh in the remaining regions to efficiently resolve the important flow 
features. One-equation and two-equation turbulence models are incorporated into 
the present unstructured approach. Results are presented for a wide range of flow 
problems including two-dimensional multi-element airfoils, two-dimensional cas- 
cades, and quasi-three-dimensional cascades. This approach is shown to gain flow 
resolution in the refined regions while achieving a great reduction in the computa- 
tional effort and storage requirements since solution points are not wasted in regions 
where they are not required. 
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1. INTRODUCTION 


The flow in axial and radial compressors and turbines is complex both in terms 
of the geometry of advanced designs and in terms of the flow structures that are en- 
countered. New designs incorporating radically shaped blades, splitter plates, and 
even multi-element blade configurations lead to even more complex flow problems. 
Early analysis of these flow problems was done mainly by experimental methods. 
Modern experimental facilities utilizing advanced measurement techniques can pro- 
vide detailed flow information for these flow problems. However, it is costly to 
perform experimental studies in the design process. On the other hand computer 
technology has seen a rapid development and the cost has decreased dramatically 
over the last two decades. This has led to significant developments in the area of 
computational fluid dynamics (CFD). A variety of solution procedures have been 
proposed for the solution of the Euler and Navier-Stokes equations. Many of these 
methods are restricted to relatively simple geometries and are difficult to employ in 
turbomachinery applications. It is therefore of prime interest to establish a general 

solution scheme for turbomachinery applications. 

The present research develops a general solution adaptive method for geomet- 
rically complex domains and complex flow structures. This approach provides a 
flexible framework for turbomachinery applications. It is important that this ap- 
proach is accurate, efficient, and easily applicable to a wide range of designs. A 
general solution adaptive method involves a combination of mesh generation tech- 
niques, solution algorithms, and solution adaptive techniques. In the following 
section recent developments in these key areas are reviewed. 
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1.1 Background 

In the area of CFD a majority of the flow solvers have been developed for body- 
fitted structured meshes. Efficient algorithms can be achieved using the body- 
fitted mesh line information. Many fast and efficient solution procedures have 
been proposed to the solution of flow for Euler [2, 14, 33, 37, 38, 65] and Navier- 
Stokes [19, 23, 53, 92] equations on structured meshes. However, applying these 
schemes to turbomachinery applications is difficult due to the problem of generating 
a structured mesh within a complex domain. The problem lies in the generation of 
a global body-fitted mesh which maps to a logical rectangle in computational space 
while satisfying a complex set of conflicting constraints in physical space. 

Two approaches have been proposed to alleviate this difficulty. One is to keep 
the structured solver and simplify the mesh generation problem. Proposed tech- 
niques include the use of Cartesian meshes [11, 12], overlaid or composite meshes [5], 
and patched meshes [68, 72]. The use of a Cartesian mesh simplifies the problem 
of mesh generation by abandoning the requirement that mesh boundaries conform 
to body surfaces. This however increases the complexity of boundary condition 
formulations in the flow solver. It can also lead to clustering of mesh in uniform 
flow regions. In the overlapping approach several subdomain grids are overlaid 
together reducing the problem to a simpler mesh generation problem within each 
subdomain. The necessary interpolation between overlaid meshes requires a special 
data structure and increases the computing time. It is also difficult to automate 
such a procedure. In a patched approach the flow field is subdivided into a se- 
ries of simpler subdomains with mesh generation performed on each block. This 
simplifies the mesh generation problem for particular geometries, but it does not 
eliminate the problem and is difficult to automate. The final mesh depends on the 
user’s experience and skill. The work of planning block subdivisions becomes a 
major obstacle in the development of computational design tools. Although these 
techniques have been successfully used in many applications, for the complex flow 
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structures encountered in advanced turbomachinery they are difficult to generalise 

and add to the complexity of the Bow solver. 

A second approach is to totally abandon the structured mesh and to formulate 
the problem using an unstructured triangular mesh. Since any distribution of points 
can be meshed with a triangular mesh, this approach eliminates any conflict be- 
tween constraints imposed to implement boundary conditions and those required to 
resolve the complex flow structures. Two approaches are most commonly used for 
generation of unstructured meshes: Delaunay triangulation methods [31, 58, 87, 88] 
and advancing front methods [48, 64, 67]. Delaunay triangulation begins with a few 
super-large triangles covering the domain of interest. Mesh points are then added, 
one by one, with a retriangulation of the mesh. Retriangulation is performed using 
the criterion that any point can not fall inside of a circle determined from any 
other three points of the existing mesh. This particular method results in an “op- 
timal” triangulation. However, the work to perform the mesh generation requires 
OiN 2 ) operations [87] because the sorting process for the triangulation is usually 
performed over all points. In addition, this method does not provide control over 
the mesh point distribution. 

In the advancing front method, the mesh point locations are determined as part 
of the mesh generation process. Starting from an initial front defined by boundary 
segments, new points are added and triangulated into the front. This process is 
repeated until the complete domain is triangulated. Mesh parameters, such as 
mesh size, aspect ratio and stretched direction, can be specified over the domain 
in this approach. Triangles are generated using the criterion that a new triangle 
should not cross any given face of the front. This method provides a great deal 
of control over the resulting mesh distribution. For turbulent flow calculations a 
highly stretched mesh is used to achieve computing efficiency. Recently, Hassan 
et d. [30] noted that the advancing front method can only produce a maximum 
allowable mesh stretching of about 10 in order to preserve mesh quality. 
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Initial algorithm work for solving flow problems on unstructured meshes was 
done in the finite element community. In the finite element approach, numerical 
schemes are formulated on an element basis, with the element formulation generally 
requiring no information from other elements. Angrand et al. [3] and Morgan et 
al. [51] have demonstrated the use of unstructured finite-element flow solvers for 
two-dimensional flow problems. These schemes use a finite difference scheme in the 
temporal discretization and a finite element formulation in the spatial discretization. 

A second approach to solving these equations* finite-volume methods, is based 
on a discrete approximation of the integral form of the governing equations. Jame- 
son and Mavnplis [36] have demonstrated the extension of Jameson’s multi-grid 
Runge-Kutta scheme [33] on regular triangles over a airfoil. Both finite-element 
and finite- volume approaches have been successfully demonstrated for the solution 
of the Euler equations using unstructured triangles in two dimensions [58, 64, 8, 90] 
and tetrahedras in three dimensions [35, 49, 66, 78]. Extension of unstructured 
schemes to the Navier-Stokes equations has recently been done for two-dimensional 
turbulent flow problems [7, 56, 59]. Finite- volume methods solve the physical con- 
servation laws directly. For turbomachinery applications the accurate prediction of 
the mass flow is critical to obtaining an accurate solution. Therefore, a finite- volume 
formulation is more appropriate. Although the unstructured approach provides a 
simple and flexible framework for solving complex flow problems, it is computa- 
tionally inefficient because of the need for the mesh connectivity information and 
they are very difficult to vectorize. It is also very difficult to implement solution ac- 
celeration methods because these techniques often assume a structured connection 
between mesh points. 

Viscous flow problems are difficult to solve on an unstructured mesh due to 
the turbulence models by which closure of Favre or Reynolds averaged Navier- 
Stokes equations is achieved. For structured mesh solvers, algebraic or zero equation 
models are the simpliest and easiest models to implement. Algebraic models do not 
have numerical stability problems and work well for a wide range of engineering 
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applications. Unfortunately, algebraic models require length scale information to 
compute turbulence quantities. The lack of body-fitted mesh line information on 
unstructured meshes makes it difficult to implement algebraic turbulence models. 
An overlaid mesh technique has been proposed to overcome this difficulty [56, 71]. 
In' this approach local structured meshes are overlaid with a global unstructured 
mesh. The algebraic turbulence model is then solved on the structured mesh. The 
necessary interpolation between overlaid meshes requires a special data structure 
and increases the computing time. Moreover the use of local structured meshes 
restricts the flexibility of the method. In order to remove the structured mesh 
dependence, more complicated turbulence models which solve one or more transport 
equations for turbulence quantities can be used. Recently, Barth [7] has successfully 
demonstrated the Bladwin-Barth one-equation turbulence model [6] on a wide range 
of turbulent flow problems using the unstructured mesh approach. Among two- 
equation models, Chien’s low Reynolds k - c turbulence model [17] has been widely 
used in engineering applications [42]. 

Without a priori knowledge of the flow structure, neither the structured ap- 
proach nor the unstructured approach can accurately and efficiently resolve the 
flow. While in principle a global fine mesh can be used to accurately resolve any 
flow structure, such an approach is impractical and computationally expensive. 
This has led to the development of solution adaptive methods. Solution adaptive 
methods can be divided into three general approaches: mesh refinement or enrich- 
ment, mesh movement, and mesh regeneration. Each type has advantages and 
disadvantages associated with it. 

In the mesh refinement approach mesh points are added or removed from the 
solution domain either by subdivision or absorption of mesh elements. Dannenhof- 
fer [21] has demonstrated this approach for a series of airfoil problems using quadri- 
lateral unstructured meshes. Starting with an initial structured mesh, irregularly- 
shaped embedded mesh regions are generated by subdividing the cells in high gradi- 
ent flow regions. This approach has been very successful in resolving complex flow 
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structures. Although it is not necessary to have a good initial mesh, the use of a 
structured initial mesh constrains the problem. If skewed cells appear in the initial 
mesh, such properties will remain in the refined mesh. Lohner [47] on the other 
hand used mesh refinement on unstructured triangular meshes. This approach lo- 
cally enrichs the mesh by subdividing triangular mesh elements. After refinement, 
any badly-formed cells are removed to improved the resulting mesh. The mesh 
refinement approach is very efficient, but it has the disadvantage of the significant 
bookkeeping involved in keeping track of modifications to the mesh. In addition, 
in both of the above formulations the adapted mesh has discontinuous variations 
in cell length scale since subdivisions are integer divisions of the original mesh. 

In the solution adaptation method based on mesh movement, the mesh point 
connectivity is fixed and the points are moved as the solution evolves. This has 
the advantage that any existing solver can be applied with minimal modification. 
To move mesh points in the structured mesh Gnoffo [29] uses an equivalent spring 
analogy, in which the mesh edges are replaced by springs with a stiffness based on 
the local gradient of some flow property. Lohner [50] and Batina [9] extended this 
approach to unstructured meshes. The disadvantage of this technique is that the 
final mesh depends on the initial mesh connectivity. 

In the mesh regeneration method, the mesh is regenerated periodically as the 
solution evolves. This may be expensive due to work required to generate the mesh 
but has many advantages. For structured meshes, mesh points are redistributed 
using structured mesh generation techniques as the solution evolves [22, 82]. In 
practice, the number of mesh points may be fixed so it has the same advantages 
as the mesh movement technique. This technique provides smooth distribution of 
mesh lines, but it is difficult to use on complex geometries. In the unstructured 
approach, Peraire et ad. [67] introduced a remeshing process in which an unstruc- 
tured mesh is regenerated using mesh parameters determined from the most recent 
solution. This approach provides smooth variation in mesh length scale and al- 
lows dense points to be placed in high gradient flow regions. Mavriplis [57] and 
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Holmes [31] demonstrated the use of Delaunay triangulation with mesh point inser- 
tion. Since Delaunay triangulation needs to search a retriangulation region when 
a new point is inserted into domain, for efficiency it requires special data structure 

to perform local retriangulation. 

- The preceding review shows that mesh generation for complex domains can 
be easily obtained through the use of unstructured meshes. An accurate predic- 
tion of complex flow problems can be effectively achieved using solution adaptive 
methods. Even for problems which can be solved using structured meshes (e.g., 
patched or overlaid meshes), there is often a significant reduction in the human 
effort required using unstructured meshes. An additional advantage of the unstruc- 
tured approach is that it provides a convenient framework for implementing solution 
adaptive methods. This leads to a great reduction in the computational effort and 
storage requirements since solution points are not wasted in regions where they are 

not required. 

1.2 Present Approach 

The solution adaptive approach used in the present work is based on mesh regen- 
eration, where the mesh is periodically regenerated as the solution evolves. While 
this is a more computationally intensive approach, it also has many advantages. 
There is very little bookkeeping required since the mesh structure is not being 
modified. Regeneration of the mesh results in a smoothly varying distribution of 
mesh points, which in turn should give better numerical solutions. Remeshing al- 
lows the opportunity to align the mesh with flow structures which, in turn, makes 
it possible to use different mesh scalings tangent and normal to a given flow struc- 
ture. A shock wave is a good example of such a flow structure, since to accurately 
capture a shock wave the mesh scale normal to the shock wave must be small. The 
ability to align the mesh will also be very important in the extension of the present 
approach to viscous flow where the mesh may also be aligned with viscous shear 
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layers and in the incorporation of flux splitting for improved resolution of shock 
waves. 

For turbulent flow calculations the level of complexity of the models determines 
computing expense, so the turbulence models used are the one-equation and two- 
equation models. Incorporating such models in an unstructured approach is still 
a very new topic and only a few approaches have been proposed [7, 59]. In those 
proposed approaches, turbulence models are solved using an implicit scheme which 
usually requires an inversion of a large matrix system of equations. In the present 
study, the one-equation and two-equation turbulence models are discretized using 
the same explicit scheme employed on the mean flow equations. 

In summary, the key elements in the present approach are the unstructured 
flow solver, the mesh generation scheme, and the adaptive remesh algorithm with 
associated refinement criteria. Jameson’s four-stage Runge-Kutta cell-vertex finite- 
volume time-marching scheme [37] is used to solve the quasi-three-dimensional 
Favre averaged Navier-Stokes equations and turbulence transport equations. The 
convergence rate is enhanced through the use of local time stepping. The quasi- 
three-dimensional equations are chosen here because they provide a better approx- 
imation to the three-dimensional flow while retaining a two-dimensional form. In 
addition, quasi-three-dimensional equations can be simplified into standard two- 
dimensional equations. This provides a more universal and convenient model for 
general flow problems. The mesh generation method with which unstructured trian- 
gular meshes are generated is the advancing front scheme first formulated by Peraire 
et ai [67]. This particular approach has the advantage of being computationally 
efficient and also provides a convenient way of adapting the mesh distribution to 
the flow solution. For viscous flows in the present work local structured triangular 
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meshes are generated around bodies to relieve the stretching limit on the unstruc- 
tured mesh generator. 

The quasi-three-dimensional flow model, the governing equations, and the Baldwm- 
Barth [6] one-equation and Chien’s [17] two-equation turbulence models are intro- 
duced in Chapter 2. The unstructured flow solver using a multistage Runge-Kutta 
finite- volume time-marching scheme is described in Chapter 3. The stability crite- 
ria, local time-stepping, and implicit residual averaging are developed. A modified 
version of artificial dissipation for highly stretched meshes and boundary conditions 
for two-dimensional airfoil and quasi-three-dimensional cascade flow problems are 
also discussed. The mesh generation procedure for both inviscid and viscous flow 
problems is presented in Chapter 4. The solution remeshing scheme and mesh 
adaptation criteria are described in detail in Chapter 5. In Chapter 6 the numeri- 
cal results of inviscid flow problems are presented. 

Unstructured solution adaptive results for the two-dimensional Euler equations 
are presented for a model multi-element airfoil, a Sanz’s supercritical compres- 
sor blade, and a Sanz’s turbine blade. Computed solutions are compared to the 
analytic solutions. Quasi-three-dimensional Euler solutions are illustrated for the 
NACA Rotor 67 transonic fan operating at peak efficiency and the Allison tandem 
blade cascade. In Chapter 7 unstructured solution adaptive results for turbulent 
flow problems are presented. Two-dimensional turbulent flow solutions for a sub- 
sonic flat plate, a RAE2822 airfoil, a NLR two-element airfoil, and a VKI turbine 
blade configurations are presented. Numerical results are compared to the available 
experimental data. A quasi-three-dimensional Favre averaged solution is presented 
for the Allison tandem blade cascade. 



2. GOVERNING EQUATIONS 


This chapter reviews the quasi-three-dimensional flow approximation for turbo- 
machinery applications. The governing equations are derived from the conserva- 
tion of mass, momentum, and energy for viscous flows on a blade- to- blade stream 
surface. The mean flow equations and turbulence models for turbulent flow cal- 
culations are also presented. With the use of the equation of state, the constant 
Prandtl number approximation, and viscosity models, the complete set of govern- 
ing equations is obtained. Initial and physical boundary conditions for the flows of 
interest are also described. 

2.1 Quasi-Three-Dimensional Flow Model 

In general, the flows in axial, radial, and mixed-flow turbomachinery designs 
are highly three-dimensional. In order to solve these three-dimensional flows in a 
relative simple manner, Wu [91] proposed the following simplification which allows 
these types of flows to be analyzed two dimensionally, but with more information 
provided. In this model, the three-dimensional flow field is described by the combi- 
nation of two separate two-dimensional flow fields as sketched in Figure 2.1. These 
two separate flow fields are composed of surfaces located in the blade-to-blade di- 
rection (Si surfaces) and surfaces lying in the hub-to-tip direction (S2 surfaces). 
In practice the flow is further assumed to follow an axisymmetric streamsurface 
as shown in Figure 2.2. The solution to the coupled two-dimensional flow fields 
requires iteration between the solution of the throughflow problem on a mean S2 
streamsurface and several quasi-three-dimensional solutions on axisymmetric Si 
streamsurfaces, since the solution for either surface requires the knowledge of the 
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shape of the other surface. lu the present research, only the flows within the quasi- 
three-dimensional SI surfaces are investigated. The radial location and thickness 
of a streamsurface is assumed to be known as a function of streamwise distance 
along the surface. This information is obtained from an ajdsymmetric through-flow 

analysis (see [41] ). 

2.2 Quasi-Three-Dimensional Navier-Stokes Equations 

The quasi-three-dimensional viscous compressible flow along the SI streamsur- 
face is expressed in terms of an axisymmetric coordinate system (m, 0) (see Fig. 2.2) 
which rotates with the blade row and is given by 

dm 2 = dz 2 + dr 2 C 2 - 1 ) 

and 

6 = 0* -Sit ( 2 - 2 ) 

where V is fixed in space and fl is the angular velocity of the blade row. In this 
coordinate system, with radius r and streamsurface thickness h taken as known 
functions of m, the Navier-Stokes equations may be expressed in the following 

form [18]: 

dtrhO) d(rkF)^d(h3) MrhS) d(h§) ( 2 . 3 ) 

- L & _ + “aS'' + _ a« 1 am + at ’ 

where 





Figure 2.2 Quasi-three-dimensional coordinate system (m, 6 ) and streamsurface 

(Katsanis [41] and Wang [86]). 
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where 

Wg = Vg- rft 

\ dr \t 1 ^ \ 

K 2 = {pVt + P - *22X7^) + ^ " ff33 *h dm 

. . „ „ E V and Vg are density, pressure, energy, m- and d- ab- 

In these equations p, p, Vm, an 9 ... 

. resnectively and W* is the relative tangential velocity, 
solute velocity components, respectively, 

The pressure p is defined by the equation of state for a perfect gas 

p = (7 - iy[£ - + 'f )1 


(2.4) 


where 7 is the ratio of specific heats. The viscous terms in the energy equal, on are 
defined as follows: 


iL, = (K r d m T + V m <7 n + Vg<7 U ) 
S 4 = {—dgT + V m <Tl 2 + Vg(Tn) 


(2.5) 


The coefficient of thermal conductivity, may be expressed as: 


C r p 

P r 


The shear stress terms are given by: 


a u = 2pd m V m + A V 'V 

<r B = 2p(»,V» + V.^)/r + * V V 

-33 = **.(!£) + 

\ *ru\ 11 are the two coefficients of viscosity which based 

In the above expressions, A and pare the two 

. v 9 /0 rp ; s the static temperature, C p is the spec 
on Stokes’ hypothesis, A = -2p/3, 1 is tne sta f 

, • Prandtl number. The dilitation is given by 

heat at constant pressure, and p T is the Prandtl numoer. 


- — 2p,o v v A±_ + lii) + - dgVg 1 

A V V ~ "*■ Vm Vdm + h dm r 


( 2 . 6 ) 
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The molecular viscosity is obtained using Sutherland’s viscosity law [76]: 

CiTl 


P = 


T + C 2 


where the constants p T , C x and C 2 for air at moderate temperature are: 
p r = 0.72, C x = 1.458xl0 _6 kg/(m - s - K 1 ' 2 ), C 2 = 110.4K 

The governing equations may be nondimensionalized with respect to chosen ref- 
erence quantities. In this work a reference length, L, and reference flow properties, 
Poo, Ko, Toe, and ^oo, are used to define the nondimensional parameters. The 
normalized variables will be expressed here as bar variables, which are 


t = 


j - m _ r - h 

L/Voo' m ~L' r = I’ h = J 


5 = — V — — V — V* n 
P Vm Kp* rn = TT 

p- - P f m JL e 

P° T ~ C PooVg' ^ 


oo 

p 


The nondimensional equations resulting from substituting the above nondimen- 
sional variables into Eq.(2.3) are similar to their original dimensional form except 

a constant, Re Lt appears in the viscous terms. Dropping the bar notation, we may 
rewrite the governing equations as: 

gfrhg) d(rkP) d(hG) ,9(rhR) d(h3), . 

at + am + sr ~ **<■ (-sr + -kr) = rhK 

where the reference Reynolds number, Re u is defined as: 

Poo V^L 


(2.7) 


Rzl = 


Poo 


and the coefficient of thermal conductivity, k t , is given as: 

P 


K-r = 


Ml - 1)AQ, 

Also note that the source term, K 2 , now becomes 

. 1 dr 


Kl ~ {pV > +p ~ + (P ~ M'vnX—) 
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For quasi-three-dimensional inviscid flow problems the governing equations are 
obtained by dropping the viscous terms in the Navier-Stokes equations. These 
equations reduce to the standard two-dimensional Navier-Stokes equations in con- 
servation form by setting r and h equal to constants. In the present work, the 
two-dimensional flow problems are solved using these equations with r = h = l. 

2.3 Favre- Averaged Navier-Stokes Equations 

The full Navier-Stokes equations provide the "exact” transport equations for 
compressible turbulent flows. With specification of appropriate initial and bound- 
ary conditions the equations may be solved directly. Unfortunately, turbulent flows 
always contain fluctuations at a wide range of frequencies and in three-dimensional 
applications it requires 0(Re»J 4 ) grid points [45] to resolve all the flow scales. Even 
with modern supercomputers direct numerical simulation of the full Navier-Stokes 
equations is still restricted to low Reynolds flows and to very simple boundary 
conditions. In most engineering applications, the mean flow properties are of pri- 
mary interest. As a result, a modified form of the Navier-Stokes equations derived 
through averaging techniques is adopted for engineering calculations. 

In the present work, the Favre- averaged technique is used to obtain the mean 
flow equations. The Favre averaging is a hybrid averaged method which uses den- 
sity weighted averaging on all fluid properties except pressure and density. For 
compressible flows this averaging results in a simpler form of the mean flow equa- 
tions than the Reynolds averaging. Performing the Favre averaging the quasi-three- 
dimensional Navier-Stokes equations (see appendix B) are 

d(rh(U )) d(rh{F)) 8 {h{G)) _ n -i Mrh{R)) + W}) ) = rk( g) ( 2 .g) 

' dt + dm .86 L dm 86 


' « 1 



MK, 

, (F) = 

wkH+m 

(pW, 

(p)V m V 9 T 

(p)E ) 


^ Vm((p)E+(p)) ) 


where 
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(p)W 9 

(p)V m W 9 

((p)V e W 8 + (p))r 
W e ((p)E+{p)) + rCl(p) 

0 

ZTi - (pV£) 

(y- [pVW))r 

V m (ZTi - (pV*)) + V e (Z? 2 - (pV"Vf)) + Kr%- (pV"ff") 

0 

- <pv"v s ") 

(SZ-0>Vp))r 

V m (Z?2 - (pVW)) + vfe - (pVp)) + fcfj - ( P V?H") 

0 

{g) = M? + o»> - - W)1 (Ji) + (<p> - p* - m) (rt) 

0 

V 0 

The equation of state for perfect gas is 

<P> = (7 - 1)W - 1 05 + i? + iff + vf n )\ 

The Favre-averaged equations are similar in form to the full Navier-Stokes equa- 
tions except for the presence of unknown stress terms in the mean momentum and 
energy equations. Among these new terms, (pV^ 2 ), (pV£V e "), and (pVg" 2 ) are the 
so-called Reynolds stresses and (pV^H") and { pVJ'H ") are the turbulence heat flux. 
These terms must be modeled to provide closure. 
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2.4 Turbulence Models 

The Favre-averaged equations contain correlations of turbulence quantit.es. In 
order to solve for the mean flow properties these terms must be modeled. Vari- 
ous approaches for modeling the correlated terms have been reviewed in the lit- 
erature (44, 85]. Although the Reynolds stress model provides general one-point 
correlation approach for turbulence quantities, it increases the level of complexity 
of equations and stiU additional unknowns need to be modeled. In the present work 
an eddy viscosity model is adopted. Using an eddy viscosity model the averaged 
equations are identical in form to the full Navier-Stokes equations, with the viscos- 
ity and thermal conductivity replaced with the corresponding effective values, W 
and (?/j>r)e, respectively. 


2.4.1 Eddy Viscosity Hypothesis 


In the eddy viscosity hypothesis turbulence quantities are related to mean flow 
properties through the use of the Boussinesq assumption, which assumes that the 
Reynolds stress tensor is linearly proportional to the mean strain rate tensor. In 
Cartesian coordinates this gives 


(~pV,"Vj') = V 



2 m 

j’dxt 



(2.9) 


In the above expression, V is the turbulent viscosity, is the Kronecker delta 
function, and k is the mean turbulent kinetic energy. Analogous to the Boussinesq 
assumption, a linear relationship between the turbulent heat flux and the enthalpy 


gradient is also adopted: 


(pV?H' f ) = 


Ji t dH _ C p ft t dT 


( 2 . 10 ) 


Pr, dXi Pr. dXi 

where Pr t is the turbulent Prandtl number, which is assumed to be 0.91 m the 
present work. Thus, the mean flow equations for turbulent flow can be expressed 
in a form similar to the original equations of motion by replacing the viscosity 
and thermal conductivity with the corresponding effective value defined as the sum 
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of laminar and turbulent parts. The effective viscosity and conductivity may be 
written as: 

l^e /^Umimr t /^turbulent 

(/*/P>~)c = (/*/ Pt ) l A min « r ~t~ (/* /Pr ) turbulent ( 2 . 11 ) 

The remaining turbulence closure problem is then to model the turbulent viscosity, 
fit, and/or the mean turbulent kinetic energy, k. 

Eddy viscosity models employing the mixing length theory [81], which is anal- 
ogous to the molecular kinetic theory, give the turbulent viscosity, Jit, as 

fi t =C(p)ql ( 2 . 12 ) 

where q and l represent the velocity and length scales of the turbulence. Depending 
on the type and the number of equations employed to evaluate these turbulence 
scales, the modeling equations are classified as algebraic, one-, and two-equation 
models. Although algebraic models are difficult to implement using unstructured 
meshes, the Cebeci-Smith [16] algebraic model is described here for complete- 
ness. In the present work, the Baldwin-Barth [6] one-equation and Chien [17] 
low Reynolds number two-equation models are employed. 

2.4.2 Algebraic Turbulence Models 

In algebraic turbulence models the turbulent velocity and length scales are mod- 
eled algebraically. The Cebeci-Smith [16] model is one common example. In the 
Cebeci-Smith [16] model, the velocity and length scales of turbulence are deter- 
mined using a two-layer model. In the inner layer where y < y c , the required scales 
and turbulent viscosity are obtained as: 

a = 


/ = KyD 


|ft| = 


hf® - S + (-- —)2 + ^_ 5^2 

y dy dx> +[ dz d y ) +[ dx dz ] 


dz 


(2.13) 
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where JC = 0.41 is the tod Kerman constant and D is the van Driest damping 
function which is used to account for near-wall viscous effects on turbulence. 

V= \ — exp(— y + M + ), A* =26 

The product of l|fi| gives the velocity scale of turbulence. The wall distance unit, 

y + , is defined as: — 

_ Ay,aUVT wa u//?waU (2.14) 

y - I'waU 

where Ay.ni is the minimum physical distance from an interior node to the wall, 
T.m, P-m -"n «* the '* al1 sheaI stress ’ deIlsity *“'* kin<!tiC visc0sity ** the 

corresponding wall location, respectively. 

Tn the outer layer, y > the turbulence quantities are given as: 


K = C a (p)U e Si 

UJi = jT°(l - u/U e )dy ( 2 - 15 ) 

where C a = 0.0168, 6 l is the displacement thickness, and U e is the edge velocity. 

The condition which defines y c is the continuity of the turbulent viscosity. In 
practice, the turbulent viscosity is determined as the minimum of the inner and 

outer turbulent viscosities. That is 

fit = min(/2‘ t ,/it) (2-10) 


Algebraic models do not require any transport equation for turbulence prop- 
erties and therefore are the simpliest and easiest models to use. However, it has 
been noted [69] that algebraic models do not account for any transport and history 
effects of turbulence and hence are not adequate for complex separated flow prob- 
lems. Incorporating zero-equation models into a unstructured approach requires 
a locally structured mesh to provide the length scale information (e.g., <5i u» the 
Cebeci-Smith model). 
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2.4.3 One-Equation Models 

Early one equation models, including those of Bradshaw et al. [13], Rubesin [73], 
and Mitcheltree et al. [62], were based on the transport equation of the turbulent 
kinetic energy, k, from which the velocity scale is evaluated. The turbulent length 
scale is then determined using an empirical description. Although the transport 
equation of the turbulent kinetic energy provides transport and history effects of 
turbulence, a length scale specified algebraically is not adequate for general flow 
problems. In the past these models have been more difficult to code and often only 
marginally better than algebraic models [69]. Further, the need of algebraic length 
scales limits the applicability of these models to unstructured meshes. 

Recently, Baldwin and Barth [6] proposed a one-equation model which involves a 
transport equation for a turbulence field variable. Except in the near wall region this 
variable is proportional to the turbulent viscosity so algebraic length scales are not 
required for the evaluation of the turbulent viscosity. This model has been shown to 
be a significant improvement over algebraic models. Although this model accounts 
for near-wall viscous effects on turbulence it does not require a very fine mesh 
spacing for resolving the viscous layer. Baldwin and Barth noted that their model 
only requires a wall mesh spacing of y + < 3.5, which relieves the numerical stiffness 
problem. In addition, since this model does not require an algebraic specification 
of the turbulent length scales it is applicable to unstructured meshes. A complete 
derivation and discussion of this model can be found in Reference [6]. The resulting 
model is summarized below. 

The turbulent viscosity is modeled as: 

v\ — C tl 'D]T>2vRT (2.17) 

In Reference [6] the transport equation for the turbulence variable, vR T , is given 
as: 

- C t ,)\JvR T V +( i7 +^ i ) V 3 {vRt) - V • V(vRt) (2.18) 
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where “P is the production of the turbulent kinetic energy. 

(dVi dVj\ dVi 2_(dV k X 1 

The corresponding nondimensional quasi-three-dimensional form is obtained as 

^ +v j!qM +i 


dm 

d{rh(RB )) | WSb)) 


dm 


dB 


do 

1 \dv t 


dv t . 


-7, 

(2.19) 


1 


where ~ x , o/~5 \ 

/ z> \ 9 (uRt) (S ) _ I djuRr ) 

<• Rb ) = am ’ ' 5b) “ r 00 

{K B ) = {C t ^ 2 -C tl )\lvR T V 

The following functions are employed to account for near- wall viscous effects on 
turbulence. 

- = (C,, 

Pi = 1 - exp(-y + /A + ) 

p 2 = 1 -exp(-y + /Aj) 

Mv*) = + + 


+ 




' ky+ 


Dj = exp(-y + /A + ), # 2 = ^exp(-y + /A+), y + = ^ 

where the modeling constants are given as: 


T wa U Ay W aU 

(^)wall ^wall 


K = 0.41, C (1 = 1.2, C (3 = 2.0 
C M = 0.09, A + = 26, At = 10 
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2.4.4 Two-Equation Models 

In most one-equation models the length scale of turbulence is described alge- 
braically. It has been recognized that the use of an algebraic length scale is not 
adequate for general flow problems [44, 85]. In order to eliminate an algebraic de- 
scription for the turbulence length scale, an additional transport equation for the 
length scale of turbulence has been proposed. A variety of two-equation models can 
be found in review articles [44, 85]. The most popular model among them is the k-e 
model. Both high Reynolds number and low Reynolds number forms of this model 
have been used successfully in engineering applications. The high Reynolds number 
k — e model is not valid in near-wall regions, so instead of integrating the mean 
flow and k — c equations up to the wall, wall functions are employed to estimate 
near-wall flow properties. A detailed numerical implementation of this approach 
can be found in References [46, 59]. This approach removes the need of extremely 
fine meshes near the wall and has been applied to a large class of flows. It is only 
appropriate to use wall functions on flows where the logarithmic law of the wall 
region exists. 

In complex flow problems the near- wall distribution of mean flow and turbulence 
properties might be different from the logarithmic law. To accurately predict the 
near-wall viscous effects on turbulence, it is necessary to integrate the mean flow 
and turbulence equations to the wall. Chien [17] modified the standard k — e model 
by accounting for the near-wall viscous effects and proposed a low Reynolds number 
k-e turbulence model which is valid up to the wall. Although this model requires 
a mesh wall spacing (y + < 1), it is coarser than the spacing (y + < 0.2) that the 
Jones-Launder [39] low Reynolds number k — e model requires. Chien ’s model has 
been widely used in engineering applications. This model is outlined below. 

The nondimensional turbulent viscosity is expressed by 

k 2 

y>t 


iirrr 


( 2 . 20 ) 
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where k is the turbulent kinetic energy 
turbulent kinetic energy. 

k=* 


and £ is the isotropic dissipation rate of the 

w v rm 


(/>)’’ ■ w 

The averaged transport equations for I and £ quantities in quasi-3D compressible 
flows are described by 

afrfc(&,» awfij) g(MgO) = fe -. ( 

— ar~ + am 3« V 3m d> ,/ 


(2.21) 


where 


(fto = | ^ 1 . <*•> = 


(p)v m k 

</>}Kne 


{&<> = 


<p)W,k 


i* <*•> 

(£ + vj dm 


l (* + *)}#) 

\(t + V l r§J 


_ V-R*L(p)t + £± 

{K kt ) -yi (frCiV - ReL^Mp)*) + & 

In Reference [17], Chien proposed the following constants and functions for this 

model: 

a k = 1.0, <r t = 1-3 


2 

9 


C M = 0.09, 

*-k 

ii 

i— 1 

& 

<4 

II 

;F M = l-exp(-0.0115y + ), 

^1 = 1, ^2 = ^ 

_(dv l dVj\ 

v -'“W i + ai) 

2 — (dVkV 

dxj 3 M ‘ V^x k ) 

„ 2/xk 

k A vl* 

S t = --r~T- ex P( 
Ayi»n 


— 9WU1 

(p)k 2 . / r««n Ay wa ii 

= Re I __ , y - y (^) wliU j5 wmU 

The extra terms, f k and &, are used to account for near-wall viscous effects on 
turbulence. Because near-wall effects needs to be resolved, this model requires a 

wall spacing y + < 1- 


2.5 Initial and Boundary Conditions 


The mean flow and turbulent transport equations presented in the preceding sec- 
tions represent an initial-boundary value problem. In order to solve these equations 
it is necessary to impose initial and boundary conditions. Numerical implementa- 
tion of these conditions will be given in the next chapter. Since only steady state 
solutions are of interest in the present work, initial conditions are not relevant to 
the solutions. 

2.5.1 Physical Boundary Conditions 

In the present work solutions for 2-D airfoils, 2-D and quasi-3D cascades are 
presented. For these cases the boundaries are either inflow/outflow or solid wall 
type. On solid walls the physical boundary conditions are the same for all problems. 
For inviscid flows the solid wall boundary condition is flow tangency condition. 

(V ~ k\r«n) • n waU = 0 (2.22) 

For viscous flows the physical boundary condition is the no-slip condition with 
either specified wall temperature or heat flux. That is, 

& & &T dT. 

T = r ~“ or *r = fel- < 2 - 23 ) 

At inflow/outflow boundaries the physical boundary conditions are dependent 
on the type of problem solved. For the 2-D airfoil problem, the inflow/outflow 
boundary is the far field boundary at infinity. The corresponding physical bound- 
ary condition is the asymptotic state of a uniform freestream condition for either 
inviscid or viscous flow. For the 2-D or quasi-3D cascade case, the physical bound- 
ary conditions at the inlet are specified total pressure, total temperature and whirl 
speed rVg. The physical boundary condition for the outflow boundary is specified 
static pressure. 
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2.5.2 Turbulence Transport Equations 

When a turbulence model based on differential equations of turbulence quanti- 
ties is used, it requires additional boundary conditions for each equation. The phys- 
ical boundary conditions are determined by the type of turbulence model used. For 
the Baldwin-Barth one-equation model, the physical boundary conditions are [6]: 

• Solid Walls: Specify Rt = 0 

• Inflow ( V • n < 0): Specify Rt ~ ( Rt)oo < 1 

• Outflow (V • n > 0): Extrapolate Rt from interior values 

For Chien’s low Reynolds number k - c model, the physical boundary conditions 
required for k and e are [43]: 

• Solid Walls: Specify k = 0, and c = 0 

• Inflow (V • n < 0): Specify k = koo, and t — too 

• Outflow ( V • n > 0): Extrapolate k and c from interior values 

It is important to realize that the total dissipation rate is not negligible at walls, 
even thought the isotropic dissipation goes to zero. The extra terms in the Chien 
k — e equations can be used to account for a nonzero dissipation rate near the wall. 
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3. UNSTRUCTURED FLOW SOLVER 

This chapter presents the spatial and temporal discretizations of the govern- 
ing equations. The Jameson’s four-stage Runge-Kutta cell-vertex finite-volume 
time-marching scheme is used to solve the governing equations. The stability re- 
quirements for this scheme are reviewed along with some convergence acceleration 
techniques. Finally the boundary conditions are formulated for two-dimensional 
airfoil and quasi-three-dimensional cascade cases. Initial conditions used in the 
quasi-three-dimensional cascade cases are also discussed. 

3.1 Finite- Volume Spatial Discretization 

For turbulent flow calculations the turbulence transport equations and mean 
flow equations may be solved in either a coupled or a decoupled manner. Although 
the coupling approach gives a single system of equations which simultaneously 
governs the mean flow and the turbulence properties, it is noted in Reference [42] 
that there is no advantage to numerically coupling these equations in terms of 
either convergence rate or accuracy. On the other hand, the decoupled approach 
allows the two set of equations to be treated separately using different numerical 
schemes. In addition, different turbulence models can be easily incorporated into 
the mean flow solver without major program modification. Since the objective 
here is to develop a flexible and robust solution scheme for complex flow problems, 
the decoupled approach is employed. With the decoupled approach the mean flow 
equations are marched in time using turbulence quantities frozen from last time 
step. The turbulence transport equations are then integrated in time using frozen 
mean flow properties. 


nr i 


** X 


Figure 3.1 Control volume for the cell- vertex scheme. 

Both the unsteady quasi-three-dimensional mean flow equations and turbulence 
transport equations are solved on an unstructured triangular mesh using Jameson’s 
four-stage Runge-Kutta finite- volume time-marching scheme [58]. In the present 
work, the cell-vertex finite-volume spatial discretization for triangular meshes [36, 
55] is used. The truncation error of this formulation has been shown to be second 
order for smooth grids and first order for general irregular meshes [70]. 

3.1.1 Mean Flow Equations 

The quasi-three-dimensional mean flow equations can be expressed in an integral 
form for control volume V with a boundary surface dV as follows: 

it JJL^+L^' 0))-nds = yy^(#)dv (3.1) 

In the cell- vertex finite- volume formulation flow properties are stored at the mesh 
points (triangle vertices) with the control volume for each point t being the union of 
all triangles with a common vertex at point i as shown in Figure 3.1. The volume 
integral is then approximated using the mean value theorem, and the boundary 
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integral is approximated using trapezoidal integration over the bounding surface. 
This process results in the follow system of semi-discrete equations for each point t 

^ + *((£?.)) = 0 (3.2) 

where R((Ui)) is the residual 

mb)) = - (tf) (3.3) 

and 

Ni Ni 

V, = '52(rhm) k A0k = - ^2(rh6) k Am k 
1=1 *=1 

Q({Vi)) = E [( r M^))*A0 t - (A(G)) 4 AmJ 

fc=l 

- fel l E [(rA(H))»Atf 4 - (A(S»*AmJ (3.4) 

k=l 

In the above expressions the subscript k refers to the k th segment of line bounding 
the control volume, Ni denotes the total number of vertices of the control volume, 
is the area of the control volume, A 0 k and Am k are the increments along segment 
k, and the components of flux vectors, ( F k ), (5*), (Rk) and (Sk), are taken as the 
average of the nodal values for each end of the segment. An entire flux balance 
throughout the flow field can be computed in one simple pass over all edges using 
an edge-based data structure. The velocity gradient in (R k ) and ( S k ) are computed 
using Gauss’s theorem [61]: 

X = iv ^ dS ( 3 - 5 ) 

where n is the unit outward normal vector to dS. Using the same control volume 
V, the above equation can be expressed as: 

^ = = (3.6) 

Replacing <j> with V m and Vg gives the velocity gradient at each node. This completes 
the discretization for the mean flow equations. 


ii* 'i 
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3-1-2 Turbulence Transport Equations 

Differential transport equations are integrated in time using the same explicit 
finite volume Runge-Kutta time integration scheme as the mean flow equations. 
This is accomplished by first taking a volume integral of the transport equations. 
Where possible the divergence theorem is then used to transform volume integrals 
into corresponding boundary integrals. The resulting integral equations are then 
discretized in space using the same cell-vertex finite-volume approximation as used 
for the mean flow equations. 


3. 1.2.1 Baldwin-Barth One- Equation Model 

The finite-volume discretization of the Baldwin-Barth one-equation model, Eq. (2.19), 
involves two gradient operators, V(IR T ) and V*. results in a complicated 

discretization. In order to simplify the discretization for the one-equation model, 
the last two terms on the right-hand side of the equation are reformulated using 

vector identities. This yields 

(v + £) V* {my - - v * ■ • VP®) = ( 5 + 2 £) V 2 (SK?) - ^ V * V (Mr) 

V oj (3.7) 

Substituting Eq. (3.7) into Eq. (2.19) and integrating over volume V gives the 

following integral form: 

v h. v {VRi) - ftdS -(* + 2 ^) / 9V v(w£) • * dS 

= (3 ' 8) 

Using the mean value theorem for the volume integrals and trapezoidal integration 
for the boundary integrals yields 

+ R{v%)i = 0 ( 3 - 9 ) 

dt 


where the residual, R{uRt) « I s expressed as: 
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R{vR T )i = (V m (R B ) + W e (S B )). - (C e ,f 2 - C ei ) yj C^Vi V 2 S( vR T )i 

-nkvi { ( 5 + 2 S), S <{Rb) ThAe ~ (Sb) kAm)k 

l ffi ) 

({ R B )rhA9 - (S B )hAm) k | (3.10) 


vt — Cp'DiDziyRT') 


In the above expression the production term, T* , is simplified using an approxima- 
tion suggested in Reference [26]. In Cartesian coordinates, it is described as 

_ ~2 ~ ,du dv.~ 

v = " 5 = v '% + ai> ( 3 -n) 

In the quasi-three-dimensional coordinate system (m, 0), S is given as 




l dV m _ (1±\ 

r 86 9 \r dm ) 


(3.12) 


The values of (R B ) and ( S B ) are computed by replacing <f> with (VR T ) in Eq. (3.6). 
For the Baldwin- Barth one-equation turbulence model, wall units or wall dis- 


tance is used to account for near-wall viscous effects on turbulence. In addition 
to employing mean flow quantities, the wall distance, Ay waa , or wall units, y+, on 
which damping functions depend have to be evaluated in order to compute the 
turbulence quantities. 

The definition of the wall units is given by 


_i Aj/w«Hv/ T wail/Pwail 

y = 1 

^wall 


(3.13) 


where Ay waU is the minimum distance from a vertex point to the wall. Other 
properties are evaluated at the nearest wall location. Based upon the procedure 
suggested by Barth [7], the evaluation of y + is implemented as follow: 


Step 1. Compute the minimum distance from each vertex point to boundary edges. 

Step 2. Store information concerning the corresponding boundary edge and weight 
factor for interpolation of the physical quantities on the boundary edge. 


nin'-r r- 
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The local wall shear stress, density and viscosity can be interpolated using informa- 
tion stored in step 2. From this information and the minimum distance obtamed 
in step 1, y + can be evaluated at each point. 

3. 1.2.2 Chien Low Reynolds Number k — e Model: 

In a similar manner, the finite-volume discretization of Chien’s k - e equations, 
Eqs. (2.21), gives 

+ R{(Cu))i = 0 ( 3 - 14 ) 

dt 

where R((Uke))i is the residual of the k - c equations 

= QUMk _ Rel'(K kl ), (3-15) 

and 

Q((U k ,))i = i;[( r Ma,))sA«a-(/.{^))aAms] 

k=\ 

- Re~L E [(rh{Rkt))kMk - (MSk<»*Am fc ] (3.16) 

fc=l 

The definitions of {&*)* <&«>*, <&>*» ^ <*“>’ ° “ be 

in Eqs. (2.21). The wall distance, Ay w «n, and wa ^ units, y + , which are used to 

account for the near-wall viscous effect on turbulence are evaluated using the same 

procedures described for the Baldwin-Barth one-equation model. 

3.2 Artificial Dissipation 

The finite-volume spatial discretization used here is equivalent to a central dif- 
ference scheme. As with other central difference schemes applied to inviscid flow 
problems, additional artificial dissipation is required to damp out the high fre- 
quency error modes. In viscous flow computations, artificial dissipation is required 
in convective dominated regions of the flow field. Even though the viscous terms 
provide physical dissipation in the viscous flow regions, a small mesh length scale is 
required in all directions to provide enough resolution for the shear stress terms. In 
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practical applications fine mesh spacing is only employed in the direction across to 
the boundary layer, so the artificial dissipation is still necessary in the str eam wise 
direction for viscous flows. 

In inviscid flow calculations, the adaptive artificial dissipation introduced by 
Jameson and Mavnplis [36, 58] is used to suppress odd-even point decoupling. 
This model is composed of a blend harmonic(Laplacian) and biharmonic operators 
which provide a weak but sufficient dissip ati ve term throughout the smooth regions 
of the flow and a stronger dissipation to suppress oscillations in the regions near 
shock waves. 

For the control volume shown in Figure 3.1 the conservative form of the dissi- 
pation operator is given by 


mm 

V 2 0i) 


23 Aj, - { Ui )) + «J(v 2 {^i) - v 2 (#))l (3.17) 

j= 1 


Ni 




1 


= Em,) - m) 

;=i 

= (Aj, + A,)/2, A,- = 


Vi_ 

At? 


_ (to) - Os)) 


= € 


i Jmm + o*» 


, £j = max(0, e" - c)) 


In the above dissipation formula, A,,, the average spectral radius of the Euler Jaco- 
bian matrix at the cell face is estimated using the value of V/Af*, V,- is the volume 


of the control volume, A t is the time step limit for a Courant number of unity, and 
d and d are constants. A Jt is chosen to provide a proper scaling in these equations. 
e j is proportional to a Laplacian of the pressure, ej and ej provide the desired 
switch between the harmonic and biharmonic operators. 


In viscous flow computations, high aspect ratio meshes are usually used to re- 
solve shear layers. Because the mesh length scales in the normal and tangential 
directions to the body are so different, the nearly isotropic scaling, A j,-, used in 
inviscid flow calculations may produce excessive dissipation in the tangential di- 
rection and therefore reduce accuracy of the viscous calculations. To alleviate this 


HI T — 
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problem, an anisotropic scaling has been introduced to provide different scales for 
different mesh coordinate directions on high aspect ratio cells. 


3.2.1 Eigenvalue Scaling 

For structured mesh calculations, the eigenvalue scaling approach [15] gives 


A* = < 


A € , for ( - direction 
A„, for Tj - direction 


(3.18) 


where A* and A, are the maximum eigenvalues for the Euler Jacobian matrices in 
each mesh coordinate direction: 


A € = (|V<| + a)A^ A„ = (|%| + «)A* (3.19) 

Since there is no apparent mesh coordinate direction for unstructured meshes, a 
scaling is performed on each edge of a control volume in order to extend this ap- 
proach to the present unstructured scheme. Consider a local coordinate system 
(n,s) on which s is tangential to the edge (ji), and n is normal to the edge. Then 
the directional eigenvalue along edge (ji) direction can be estimated as: 


A. = (|V.| + a)An 


( 3 . 20 ) 


In the above expression, V, is the s-component velocity, An is the length scale in 
the tangential direction to edge (ji), and a is the sound speed. Using information 
at points i and j we can estimated all the required parameters in Eq. (3.20) except 
for An. Further information is required to estimate An. 

In the edge-based data structure, an additional triangle on either side of the 
edge can be constructed, the distance between centroids of these triangles is used 
to estimate the normal length scale An. Replacing the A* with the directional 
eigenvalue A, gives the anisotropic scaling for large aspect ratio unstructured mesh 

cells. 

Although the anisotropic scaling gives appropriate dissipation on large aspect 
ratio meshes, various investigators have found that such a scaling may not produce 
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enough dissipation when mesh becomes highly stretched. This leads to slow con- 
vergence. To overcome this problem, a new scaling based on the combination of 
isotropic and anisotropic scaling was introduced by Martinelli [52]. On a structured 
mesh, the scaling becomes 

* = ** (‘ + = (‘ + (\) ) (3 - 21) 
On the present unstructured mesh, it is modified as 

A* = A, ^1 + ) (3.22) 

In the above expression, a is a constant. A n is the eigenvalue in the direction normal 
to the edge and is given as 

An = (|K| + a)As (3.23) 

where V n is the n-component velocity and As is computed as the distance between 
points i and j. 

3.2.2 Local Velocity Scaling 


For viscous computations, the physical dissipation in near wall regions produces 
enough dissipation in the direction normal to a surface, because a fine mesh spacing 
is used in this direction. This allows the artificial dissipation term in the normal 
direction to be reduced. In the present approach, a local velocity scaling is incor- 
porated into the directional eigenvalue scaling. That is 

(K 

U 

where q is the local total velocity and 0 is the angle between edge (ij) and a vector 
tangential to the nearest wall location. 

Adding the modified artificial dissipation to Eqs. (3.3), (3.9), and (3.14) results 
in the fined form of the finite-volume spatial discretization: 

mOi)) = - £((#»! - *((£/,» 


) 


(3.24) 



nr-?'!- 


(3.25) 
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The dissipation term adds a third-order error in smooth regions of the flow, and a 
first-order error near shock waves, so a second-order accurate finite-volume scheme 

is preserved except in regions close to shock waves. 

The implementation of the dissipative terms requires two loops over all the edges 
in "the domain. In the first loop a harmonic (Laplacian) operator is computed at each 
node. In the second loop the blended adaptive artificial dissipation is obtained by 
simultaneously collecting the Laplacian terms needed for the biharmonic operator 
and the terms for the harmonic operator. 

3.3 Runge-Kutta Time-Integration Scheme 

The semi-discrete equations are integrated in time using the following q-stage 
Runge-Kutta scheme: 

(U) w = { U) n 

<t/)(D = (ff)<°> - AtR{(U)W) 

: (3.26) 

(jj)ii- 1 ) = (U) {0) - 

(#)<«> = (U) {0) - a, A LR«tf) (, - 1) ) = (# > n+1 

where the superscripts in parentheses represent the particular stage of the scheme. 
(U) n denotes the value of { U ) at time step n, a u ■ • • , a, are the coefficients of the 
particular multi-stage scheme, and R{{U)) is the residual of the spatial discretiza- 
tion as given in Eq. (3.25). Runge-Kutta schemes were originally developed for 
high temporal accuracy. However, for steady-state problems, time accuracy is not 
as important as efficiency of the scheme. This has led to the use of hybrid mul- 
tistage schemes which allow the convective and dissipative terms to be evaluated 

separately. 

For the results presented here a standard four-stage Runge-Kutta scheme with 
coefficients ai = 1/4, a 2 = 1/3, a 3 = 1/2 and a 4 = 1 is used. The viscous operator 
((R), (§)) and the artificial dissipation operator £>((&)) are evaluated in the first 
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stage step and frozen for subsequent stages. In this way the CPU time required for 
the evaluation of the dissipative term is reduced by a factor of four. 

3.3.1 Stability Criteria 

. The time step for all explicit time marching schemes is restricted by some sta- 
bility limit usually expressed in terms of the CFL number. For the Euler equations 
the maximum allowable time step is described as: 


At = CFL • At c 


(3.27) 


where A t c is the convective time step. For the Navier-Stokes equations the time 
step limits [54] due to both convective and diffusive effects must be considered, 
which may be expressed as: 


At = CFL 


At c • At v 


At c + At. 


(3.28) 


where CFL is the Courant number and A t v is the diffusive time step. 

The Courant number for the multi-stage scheme is obtained by performing the 
von Neumann stability analysis for the one-dimensional model problem: 


Uf -J- u x 4- fi Ax u xxxx — 0 (3.29) 

The stability region for the hybrid multi-stage scheme depends on the discrete spa- 
tial operator, the coefficients, a,-, the number of stages, q, the manner of evaluations 
of the artificial dissipation, and the smoothing coefficient, /x. Although for a m-stage 
scheme a maximum allowable stability limit, CFL < m — 1, can be obtained through 
optimizing the coefficients, a,-, such a scheme is not necessarily optimal in terms 
of computational effort or accuracy. The present four-stage Runge-Kutta scheme, 
using the centered difference operator with the dissipation computed in the first 
stage and frozen for subsequent stages, has been shown to be stable with a Courant 
number CFL < 2.6 for the one-dimensional model problem with /x = 1 /32. Note 
that this analysis only provides a reference value of the maximum Courant number 


iirrr 
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for the present scheme. In practical applications values less than the maximum 

Courant number might be needed to assure stability. 

The values of A t e and Ai„ are obtained by performing the von Neumann sta- 
bility analysis on the finite difference equations. Since it is not possible to perform 
a stability analysis directly on the present unstructured scheme, the stability limit 
for the multi-stage scheme will be inferred from the corresponding structured for- 
mulation. The time steps A t c and &t v for a general two dimensional structured 

mesh [43] are given as: 


At c = 


Vi 


A C £ + 


A t v = 


Vi_ 

A„ t + A. 


v v 


(3.30) 


where 


A c< = (|V*| + a)Ai/, A*, = (|K,| + a)A£ 


1 

Av< " Re L Vi 


47/^ 


Av ” " Re L Vi |»(Pr)« 


(p)(Pr)e 

47/x* 


a ” j+ 4 a H 


y\ 


A£ = \[A + ^ = y/ x l + 

Although the edge-based approach provides local mesh coordinate information for 
evaluating the above equations, it is expensive to compute all the information for 
each edge. In order to reduce the work of evaluating the time steps m the present 
unstructured scheme, the above equations are simplified using time steps in the £ 
and T} directions. 

At c = min{At c< , At c „}, At* = min{At v< , At Vv ) (3.31) 

After neglecting the cross terms in A v , the time steps in the £ and tj directions 
might be described as: 

A£ „ . At/ 


Ate, = 


A t v , = RclAtj 


C( ' M + «’ 

2 (p)(P»~)e 


A ‘ e "~ |V,| + o 


47Pe 


At Vri = RclA£ 


■2 (p)(Pr)' 


47P< 
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For the cell- vertex unstructured scheme, it is not necessary to compute the time 
steps in both tangential and normal directions for every edge, because the tangential 
direction of an edge can be approximated by a normal direction of another edge 
when looping over all the edges of a control volume. In the present work, the 
time step on the direction normal to a edge is used. Therefore, the time steps 
corresponding to unit Courant number on each edge may be rewritten as: 


A '« - = 4 >( -)* Ar; n 

\V-n\+a 47 ft e 


(3.32) 


In the above expression subscript k denotes the k th edge surrounding the control 
volume, Ar n is the normal distance from interior point i to the edge, Ar„ is the 
distance between two end points of the edge, velocity V and speed of sound a are 
taken as averaged values along edge k, n is the unit vector normal to the edge, and 

• n| -+■ a)i is the maximum wave propagation speed along the direction normal 
to the edge. 

In order to advance the solution in time with the maximum allowable speed, 
the maximum CFL number is always applied in the estimation of time steps. 


3.3.2 Local Time-Stepping 


If only the steady state solution is of interest, convergence to the steady state 
solution may be accelerated by sacrificing the time accuracy and marching the equa- 
tions at each node in time by the maximum permissible time step, as determined by 
a local stability analysis. Based on the above edge-based approach the permissible 
time step at node t may be taken as the minimum of the allowable time steps from 
all surrounding edges to that node. 


A<Ci = Atvi = ™£{ A ^} 


A ti = CFLtn^x 


A tc, • Af p< 
Ate, + A t Vi 


(3.33) 


in 1 
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3.3.3 Implicit Residual Averaging 

Another simple method of increasing the allowable time step is implicit residual 
averaging. The concept of implicit residual averaging is to increase the stability 
limit of the basic time-integration scheme by implicitly smoothing the residuals. 
The implementation is performed by implicitly solving the equation: 

Ri = Ri + e V 2 Ri ( 3 ‘ 34 ) 


where R!s are the smoothed residuals, and e is a smoothing coefficient. The stability 
limit of the smoothed scheme can be estimated using the following inequality [34]: 


1 ( CFP \ 
~ 4 \CFIr 2 7 


(3.35) 


where CFIT is the Courant number of the unsmoothed scheme. 

For the unstructured formulation, a Jacobi iteration method has been sug- 
gested [58] for solving the implicit residual smoothed equation. Due to the fact 
that the computational work involved in solving Eq. (3.34) should not outweight 
the gain in convergence rate, two iterations were employed. With the same Lapla- 
cian operator as in the artificial dissipative term, the Jacobi iteration scheme gives 

™»+i _ L ± iS (3.36) 

It is important to note that the Jacobi iteration method converges slowly unless 
the matrix is strongly diagonally dominant. This implies that a small value of e is 
required and this restricts the increase of the stability limit. In the present work, 
it has been found that using two Jacobi iterations with e = 0.5 saves about 20 % 

in CPU time. 

In order to improve the performance of the residual averaging method, an alter- 
nate iterative method using an approximate factorization is employed to enhance 
the convergence in solving Eq. (3.34). In matrix form, Eq. (3.34) gives 


[Mij]Ri = Ri 


(3.37) 
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[A/ij] can be factorized using the corresponding diagonal lumped mass matrix [63], 
that is 

[Mij] = [Da] + ([Mij] - [A;]) (3.38) 


where 


N 


j = 1 


Substituting Eq. (3.38) into Eq. (3.37) gives 


[A;] SRi . Ri - [M^] R” 

A* +1 = R? + 6 FU 


The above equation can be simplified as: 

Sn+1 _ A + tNjFCj -f e R% 

1 1 + 2 Nit 


(3.39) 


(3.40) 


This method requires about the same computing work as the Jacobi iteration 
method. The advantage in this iteration scheme is that the value of t can be 
increased up to 1.0 with a resulting 40 % savings in CPU time. This method is 
used in the present work. 


3.4 Initial Conditions 


As mentioned previously, the mean flow equations and turbulence transport 
equations represent an initial-boundary-value problem. Appropriate initial and 
boundary conditions are needed to solve these equations. For steady state solutions, 
initial conditions are used only as a starting point for time-marching and only affect 
the convergence rate to the steady state solution. However, the imposed conditions 
should not be so inconsistent that the time-marching scheme diverges. 

3.4.1 Mean Flow Equations 

-For external flows initial conditions are not critical to the iteration, and uni- 
form freestream conditions are generally imposed. For turbomachinery applications, 


in 
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there is usually high loading and high turning of the flow. It is not generally pos- 
sible to start a solution with constant initial conditions. In the present work the 
quasi-one- dimensional solution introduced by Chima [18] is imposed to provide a 

smooth variation of initial flow conditions. 

' For quasi-one-dimensional flows with area change the continuity and energy 

equations which state the conservation of mass and rothalpy are 

m = prhAOW cos a = constant ( 3 * 41 ) 

/ _ c P T' - rttVe = CpT" - = constant (3.42) 

where ()' and ()" values denote absolute and relative total conditions, respectively, 
A 9 is the blade spacing, and I is the rothalpy. For flow outside the blade row the 
free vortex flow assumption gives 

rV e = constant ( 3 - 43 ) 


From Eq. (3.41) the relative total velocity can be described as 


From the isentropic relations we have 

P_ 

p" v r 


P / T . _L_ 

— = ( — V - 1 

j, V T"’ 


T w 2 

T» = 1 _ 2CpT" 

^°’ W 2 

p = p '( l ~ 2CpT"^~ l 

Using these relations Eq. (3.44) may be rewritten as 


where 




m 


p"rhA6 


( 3 . 44 ) 


(3.45) 

(3.46) 
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Within the blade row the free vortex assumption is replaced with the assumption 
that the flow angle a varies linearly through the blade row. The relative total 
velocity can be described as 

W 2 i 

w ’ - ~ 2 = 0 (3 ' 47) 

where 

<f> = — 

p"rhA0 cos a 

After the inflow conditions and the relative flow angle at the trailing edge are 
specified, the solution of the quasi-one-dimensional flow can be obtained from up- 
stream to downstream. For each point upstream of the blade the relative total 
velocity is obtained using a Newton iteration method to solve Eq. (3.44). Other 
flow properties are then obtained thru Eq. (3.43), constant rothalpy, and the isen- 
tropic relations. Within the blade row Eq. (3.47) is solved for W. Once the flow 
condition are computed at trailing edge, Eq. (3.46) is used again for the downstream 
region. 

The quasi-one-dimensional solution is not necessarily consistent with the flow 
tangency or the no-slip condition. To avoid an abrupt change of flux at a near wall 
cell, the flow tangency or no-slip condition is enforced at the blade surface and a 
smoothing operator is then used to smooth the flow quantities near blade edges. 

3.4.2 Turbulence Transport Equations 

For the Baldwin-Barth one-equation model, the turbulence field variable is ini- 
tialized using specified freestream turbulent viscosity: 

Rt = Rtoo < 1 

For the Chien low Reynolds number k — c model, the initial distribution of k and e 
are based on specified freestream turbulent viscosity and turbulence intensity. The 
turbulence intensity is defined as 



HIT I" 


(3.48) 
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In the present work, the values of %*> and v too are specified as 

%o < -0.1 %, v too < 1 ( 3 - 49 ) 

Substituting Eq. (3.49) into Eq. (3.48) yields the freestream turbulent kinetic en- 
ergy. Then, the freestream dissipation rate can be obtained as 


3.5 Boundary Conditions 

For a cell-vertex scheme flow properties are required at boundary nodes. In 
the two-dimensional or quasi-three-dimensional mean-flow equations four bound- 
ary conditions are required along each boundary of the domain. In addition to the 
physical boundary conditions mentioned in the previous chapter, additional nu- 
merical boundary conditions might be needed to close the set of equations. When 
differential turbulence transport equations are used for the closure of the Favre- 
averaged Navier-Stokes equations, additional boundary conditions are required on 
each differential transport equation. The boundary conditions depend on the type 
of turbulence models used. 

3.5.1 Inflow/Outflow Boundary Conditions 

At inflow and outflow boundaries, a nonreflecting or radiation boundary condi- 
tion based on a characteristic analysis is implemented to allow the numerical errors 
to propagate out of the domain. The idea of a characteristic formulation is try to 
correctly capture the physics of the flow in terms of a wave propagation problem by 
considering the characteristic waves in incoming and outgoing directions separately. 
For the incoming waves, the characteristic variables are held at prescribed values, 
while for outgoing waves the variables are extrapolated from the interior. 

The characteristic formulation is derived in a local coordinate system tangential 
and normal to the boundary. Assuming the normal variation to be much larger than 
the tangential variation of state variables U , the Euler equations may be expressed 
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as 


where 


dW . dW 

dt +A dn 


= 0 


A = Diag{A x , A 2 , A 3 , A*} 


' w x ' 


P - p/a 2 


’ A x ' 


9n 

w 2 




a 2 


9 n 

W 3 


(<7n + p/pa)/>/2 

5 

A 3 


9n + a 

w 4 


. (“9* + P/pa)/^ 


. A 4 


9 n - a 


(3.50) 


(3.51) 


where W, is the characteristic variable, A,- is the associated wave speed, n is the 
local unit inward normal vector to the boundary, and q n and q, are the normal and 
tangential velocities, respectively. Barred values are evaluated using values frozen 
from the last Runge-Kutta stage step. 

Along each characteristic line defined by the corresponding wave speed, Eq. (3.50) 
reduces to a set of ordinary differential equations. That is. 



along ^ = A; 


(3.52) 


Employing Eq. (3.52) along the boundaries of the computational domain gives 
characteristic formulations for solving for the flow properties at boundary nodes. 
The discrete form of Eq. (3.52) involves the characteristic variables propagating 
from outside the domain as well as from the interior. The characteristic variable 
Wi must be specified if the corresponding wave speed A j originates from outside 
the domain. Otherwise, it is extrapolated from the interior. Based on the local 
velocity normal to the boundary (see Figure 3.2), there are four possible boundary 
formulations summarized as follows: 

For subsonic inflow (0 < q n < a ), there are three incoming waves corresponding 
to A x , A 2 and A 3 , and one outgoing wave corresponding to A*. Therefore, the char- 
acteristic variables, W \ , W 2 , and W 3 are prescribed and W 4 is extrapolated from 
the interior. This leads to the following system of equations. 


n r 
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Figure 3.2 Characteristics at far field boundaries. 
Pb ~ Pb /“ 2 — Pp t ~ ft "/® 2 

9n b + Pb/P 5 = 9"pr + ft*/ P G 
-9n b +Pb/P« = -9n„+P«c/P a 


(3.53) 


where the subscripts b, pr, and ex represent the boundary, prescribed, and extras 
olated values, respectively. Solving Eq. (3.53) for the boundary values gives 


9a b 9a pr 

Pb = » [ft* + Ppr + P®(9npr — 9n««)] 

Pb = Ppr + (Pb ~ Ppr)/« 2 
9n b = 9n p , + (Ppr-Pb)/P« 


(3.54) 


For supersonic inflow (0 < a < ?»), aU four waves (A;) come inward, so all four 
characteristic variables (WJ) are prescribed. This is equivalent to prescribing aU 
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four flow properties. 


^ 3 b ^7*pr 

Pb = Ppr 

Pb ~ Ppr 

9n b = q npr (3.55) 

For subsonic outflow (—5 < q n < 0), there is only incoming wave corresponding 

to A 3 , and thus the corresponding characteristic variable W% is prescribed. Using 
the same procedure as for the subsonic inflow case the boundary values are obtained 
as 


Pb = ^ [Pe*+Ppr + M?n P r ~ <7n„)] 

Pb = Pex + (Pb Pex)/^ 2 

?n b = + (Pb ~ Pex)/pa (3.56) 

For supersonic outflow (<y n < —a), all four waves A,- propagate from the interior, 
therefore all four characteristic variables are extrapolated from the interior. 
This is equivalent to extrapolating all flow properties from the interior. 


— 


Pb — Pex 

Pb = Pex 

?n b = q^ x (3.57) 

When the prescribed and the extrapolated values are given, these relations, Eqs. (3.54) 
to (3.57), determine the primitive variables p, g n , and p at the boundary. Energy 
is obtained from the equation of state. 

For the two-dimensional airfoil case the far field boundary is, in practice, placed 
at a finite distance from the airfoil. In order to accurately approximate the asymp- 
totic state of a uniform freestream condition the far field boundary must be placed 


n i • 
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at a distance of 0(100) airfoil chord lengths from the airfoil. For subsonic flows the 
distance from the outer boundary of the domain to the airfoil can be reduced by re- 
placing the uniform freestream far field condition with the vortex far field boundary 
condition. Based on the work of Usab [83] the vortex far field boundary condition is 
defined by the combination of freestream and a compressible point vortex centered 
at the airfoil quarter chord. Use of the point vortex correction allows a reduction 

in the distance to the far field boundary by a factor of 10 [21]. 

For turbomachinery applications the total pressure, total temperature and whirl 
r y $ are specified at the inlet with exit pressure imposed at the outlet. In order 
to explicitly prescribe these quantities, an alternate boundary formulation based 
on Riemann invariants is used at the inlet boundary. At the outlet boundary a 
characteristic formulation is used. 

The Riemann invariant formulation is based on the work of Chima [18]. In 
this formulation, an upstream-running Riemann invariant is extrapolated from the 


domain interior: 

d m R~ — 


V 1 _ 


j(V? + *v«)(~) + <■’'"(^>1 ( 3 - 58 > 


where 


R- =V m - 


2a 

7 -l 


is the upstream-running Riemann invariant. 

For subsonic Bow at the inlet the Riemann invariant is extrapolated using 
backward-differencing in Eq. (3.58). Based on IT, the isentropic relation, the 
specified total temperature and whirl rV», the velocity V„ at the boundary is given 

by: ( 7 - pit- + t/( 7 TT)( 4 CpT- - 2 v?) - 2(7 - (3 59) 

v m — (7 + 1) 

Density and pressure are computed using the isentropic relations and the specified 
inlet conditions. The nondimensional specified total pressure and total temperature 
yield the total density: 

p ' = (7 MDf ; 
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Prom R~,V m , and rV$, the total velocity and speed of sound at inlet boundary are 
found: 

«b = y/V’ + V? 

«b = 2 -^-i(K,-.R“) (3.60) 

Then, the isentropic relations and equation of state give: 



'b = + f (12 + V?) (3.61) 

For supersonic inflow all four quantities are specified. 

At the outlet, the characteristic formulation with prescribed exit pressure are 
used for subsonic flow, while all four variables are extrapolated from the interior 
on supersonic flow. 

3.5.2 Solid Wall Boundary Conditions 

The cell- vertex scheme flux balance is evaluated using flow properties at node 
points. This requires boundary values to be updated at each time step. For inviscid 
flow calculations, the flow tangency condition is enforced at wall. Along solid wall 
boundaries the imposition of the no flux condition in flux balance calculations 
does not guarantee that the flow tangency condition will be satisfied. To enforce 
the flow tangency condition the solid wall boundary conditions are implemented 
in a predictor- corrector manner using a characteristic analysis. With q n = 0 in 
Eq. (3.51) the characteristic formulation shows one outgoing wave requiring one 
prescribed condition. 

The state variables are first predicted based on the Runge-Kutta scheme: 

U p = U n + AU 


'III 'TP 


(3.62) 
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Substituting the condition = 0 and the other variables from the predicted state 
into Eq. (3.51) results in a characteristic formulation of the solid wall boundary 
conditions. The boundary values are then corrected using the following formulation. 


9n c = 0 

= <1; 

Pc = Pp + ^ P pa 

p c = Pv + W/a ( 3 ' 63) 

where the subscripts p and c represent the predicted and corrected values, respec- 
tively. Again energy is computed using the equation of state. 

For viscous flow calculations, the no slip condition with either adiabatic or 
isothermal wall condition is specified at the wall. For the cell-vertex scheme an 
additional numerical boundary condition is required to update all boundary flow 
properties. In practice a zero normal pressure gradient is imposed at the wall. The 
boundary conditions imposed at the wall are 


V = V waB 
T = T wa ll, 

( ^ ) ”“ = 0 


i aT \ — o 

OT W”" 


(3.64) 


Density is computed using the isentropic relations and equation of state. 


3.5.3 Periodic Boundary Conditions 

For a single-blade-passage a spatial periodic boundary condition is imposed on 
the upper and lower boundaries between which an equal pitch spacing is maintained 
from inlet to exit. In the present approach the points are placed on boundaries 
and no imaginary cells are placed outside periodic boundaries. Along the periodic 
boundaries the flux balance is performed in the same way as at the interior points. 
The net change at a boundary point is simply the sum of the partial sums from the 
corresponding points on the upper and lower boundaries (see Figure 3.3). 
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Figure 3.3 Flux balance implementation for periodic boundary points. 

3.5.4 Special Boundary Conditions 

For turbomachinery applications special boundary conditions are required. Blades 
designed on an analytical basis provide good test cases for the validation of numer- 
ical methods. Many of these cases have an open profile at the blade trailing edge. 
Changing the blade shape to close the trailing edge often results in a very different 
solution. In order to obtain an accurate numerical solution, the nonclosed trailing 
edge needs to be correctly modeled. One example of a nonclosed profile is the "vis- 
cous” blade profile obtained by adding the boundary layer displacement thickness 
to a physical blade. The inviscid solution computed using this "viscous” profile is 
equivalent to an inviscid solution about the true blade with surface injection used 
to account for the boundary layer displacement effects. Based on this analogy the 
nonclosed trailing edge is modeled as a free jet which accounts for the boundary 
layer displacement inviscidly. 

In the present work, the nonclosed trailing edge is modeled as a uniform free 
jet at the blade trailing edge. The uniform flow properties are evaluated using the 
averaged values of the upper and lower blade surface properties at the trailing edge 
(see Figure 3.4): 


It i 


51 


Pi = (P L + P U )/ 2 

Uj = (« L + « l/ )/2 

Uj = (u L + v u )/2 

ej = (^ + 6^/2 ( 3 - 65 ) 

where the subscript j denotes the jet flow end the superscripts L and U represent 

the lower and the upper end points, respectively. Pressure is obtained from the 
isentropic relations and equation of state. 




Figure 3.4 Uniform jet assumption at the nonclosed trailing edge. 

For turbomachinery applications the exit pressure is set or measured in an av- 
eraged (mass- or area-averaged) sense. For a uniform outlet flow this is numerically 
equivalent to specifying a constant pressure at the exit. However, there are situa- 
tions in which the exit flow is not uniform. For example, a turbine blade operating 
at a high speed it is likely to have oblique shock or expansion waves which propagate 
to the exit boundary. If the axial exit velocity is subsonic, one boundary condition 
is still required. Because the exit flow is not uniform, it is physically inappropriate 
to specify a uniform exit pressure at the exit. Therefore, nonuniform exit boundary 

conditions axe required. 

In order to account for the nonuniformity the exit pressure is imposed in a 
predictor/corrector fashion. In the predictor step the state variables are extrap- 
olated from the interior and a mass- or area-averaged exit pressure is evaluated 
using the extrapolated value. The exit pressure is then corrected by subtracting 
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the difference between the specified and the predicted averaged exit pressure from 
the extrapolated value: 


— Pex it Pp»T. 

Pc = Pp + Ap (3.66) 

where the subscripts c and ex denote the predicted and corrected state, respectively, 
P*v« is a predicted averaged value, and exit represents the specified exit value. 

A variation of this boundary condition is that the exit pressure is only known at 
one boundary point. Since there is no averaged exit pressure available, the nonuni- 
fonn exit pressure is imposed using a slight modification of the above procedure. In 
Eq. (3.66), the averaged pressure difference is replaced with the pressure difference 
computed at one boundary node. 

3.5.5 Boundary Conditions for the Turbulence Transport Equations 

For the Baldwin-Barth one-equation model, the boundary conditions stated in 
Reference [6] are used. 

• Solid Walls: Specify Rt = 0 

• Inflow (V ■ fi < 0): Specify Rj = ( Rt)oo < 1 

• Outflow (V ■ n > 0): Extrapolate Rt from interior values 

For Chien’s low Reynolds number k — e model, the boundary conditions listed in 
Reference [43] are used. 

• Solid Walls: Specify k = 0, and e = 0 

• Inflow (V • n < 0): Specify k = koo, and c = €«> 

-• Outflow (V ■ n > 0): Extrapolate k and c from interior values 
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4. MESH GENERATION 


In the present work an adaptive remeshing procedure is applied to the resolution 
of complex flow problems. With a solution adaptive method based on remeshing, 
the mesh is recomputed periodically as the solution evolves. Therefore, the speed 
at which a mesh can be generated is critical to the overall performance. For this 
reason both the initial and adapted triangular meshes are generated using the 
advancing front method developed by Peraire et al. [67]. This particular scheme also 
allows a significant amount of control over local mesh properties with specification 
of the local length scale, aspect ratio and orientation of triangles. As shown in 
Reference [67] these mesh parameters provide a means of directionally adapting 
the mesh. The present unstructured mesh generation scheme follows the work 
described in References [67] and [48]. 

4.1 Advancing Front Method 

The advancing front mesh generation scheme begins with the specification of 
the domain boundaries and a background mesh on which nodal values of mesh 
parameters 6 , a and s are prescribed. Referring to Figure 4.1, these mesh pa- 
rameters define the characteristics of the elements to be generated: 6 defines the 
node spacing; s defines the stretching or aspect ratio; and a defines the direction 
of stretching. These three mesh parameters are assumed to have piecewise linear 
spatial distribution over triangular elements of the background mesh. When a new 
triangle is formed, the mesh parameters are first evaluated through interpolation 
from the background mesh. Initially, the background mesh is constructed using a 


54 



Figure 4.1 Definition of mesh parameters. 

very coarse hand triangulation of the domain (e.g. see Figure 4.2). If the domain 
is being remeshed the last mesh is used. 

The first step of the mesh generation process involves generation of mesh points 
on the boundary of the domain. This step is done in a predictor/corrector fashion 
with a spline-fitting technique providing the boundary surface location. The two- 
step process is described below with the example problem shown in Figure 4.3. 

Assume background mesh information along the boundary of interest is known 
(see Figure 4.3). The spline function is given in terms of an arc length tangent 
to the boundary ranging from 0 to r^. The coordinates (x, y) and the mesh 
parameters are expressed as function of r. Marching around the boundary, a new 
list of boundary points with spacing based on the local value of 6 is generated 
iteratively. The iteration process is described as follows. 

r i + 1 = r * d" ^V+l/2 

Starting with initial guess <5? +1/2 =. 6 it r? +1 is computed. The local value of Sf +l 
is interpolated from background mesh. Then new value of Sg*. 2 is computed as a 


and background mesh for NACA0012 airfoil 


Figure 4.2 Boundary information 


o 


: Initial distribution 


: Predictor step 
: Corrector step 



56 


mean of 6, and 6" +1 

= 0-5(<T +1 + «.) 

The above iteration is continued until convergence. In order to generate a smooth 
distribution the following restriction is imposed: 

= minCLlW,,^ 1 ) 

Then, if a local minimum 8 exists within node i and t + 1, the location r and the 
value 8 replaces the values at point i -f 1. 

The above procedure produces a smooth mesh spacing distribution if the value of 
6 is in ascending order corresponding to the marching direction. In order to obtain 
a smooth variation for a general distribution of 8, the boundary points and mesh 
parameter distribution are first predicted marching backward around the domain 
boundary using the above iterative scheme. Then, the boundary points distribution 
is corrected marching forward around the domain boundary placing points on the 
boundary using the mesh parameters obtained at the predictor step. There is no 
need for iteration in the corrector step, because the predictor step provides a good 
prediction of the boundary mesh space distribution 6. 

These boundary points define the initial front list, a set of straight line segments 
which connect consecutive boundary points as shown in Figure 4.4. Step by step 
new triangles are then added alone the front and then absorbed into the front. 
Starting with the shortest segment of the front, a new node is added within the 
domain at a point determined by the local values of the mesh parameters and 
specified mesh quality constraints (see [67]). 

{ max(<7,, 0.75<7fr) if <r, < <r fr 

(4.1) 

min(<7j, 1.3<r fr ) if <r, > a b 

where a represents the mesh parameters, 8, s, and a, and the subscripts i and fr 
denote local value and value of the active front, respectively. 

Based upon the local values of the stretching parameter, the Cartesian coordi- 
nates (x,y) are transformed into elliptic coordinates (x e ,y e ), in which the triangles 
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Figure 4.4 Initial front set up. 


satisfying the stretching conditions will look "equilateral”. This transformation is 
described as 

/«.\ r 1 /.« 0 1 T rns at sin m ]( x \ 

(4.2) 




1 js 0 


V. ) 


0 1 



cos a sin a / x \ 

— sin a cos a V y I 


where s and a are the local values of mesh stretching and mesh orientation. In 
elliptic space, a list of candidate nodes for the triangulation is created searching 
the front nodes which lie inside of a circle formed by the new point and a specified 
radius. The radius used in the present work is set equal to 1.66. The candidate 
points are then sorted in ascending order according to their distance from the new 
point. The new point is placed at the first place in the list if the distance from the 
new point to the end points of the active front are less than 1.66. Otherwise, it is 
placed at the end of the list. 

-A new triangle is formed using the active front face and the first node from 
the candidate list which satisfies a criterion that my face of the new triangle can 
not cross any existing front segment. The front list is then modified to include the 
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faces of the new triangle in the front, while removing the triangle face adjoining 
the previous front. After the new triangle is created the elliptic coordinates are 
transformed back to the Cartesian coordinates. This process is repeated until there 
are no segments in the front list, at which time the domain has been completely 
filled with triangular elements (e.g., see Figure 4.5). 

To improve the quality of the flow solutions computed on the new mesh the 
following three additional operations are performed in the present implementation. 
First each node within the domain which is a vertex for less than 5 triangles is 
removed and the mesh is retriangulated in the region of the removed node. Second, 
diagonal swapping is performed to remove any obtuse triangle. If the angle between 
any two sides of a triangle is greater than some specified angle a. w the diagonal 
of two adjacent triangles are to be swapped. In the present work a sw = 170 deg is 
used. Last, the mesh is smoothed to remove any nonuniformities in the mesh using 
a Laplacian type operator: 

X?* 1 =X? + - £(*7* - XD (4-3) 

« ft 

where X m+1 is the new location of the point Xi, n is the number of the neighboring 
points and e is the smoothing coefficient. In the present work two smoothing passes 
with e = 1 are used. 

4.2 Structured Triangular Mesh Algorithm 

The advancing front technique provides a flexible way to generate unstructured 
triangular meshes for complex geometries. However, Hassan et al. [30] have recently 
noted that the advancing front technique can only produce a maximum allowable 
mesh stretching of about 10 in order to preserve mesh quality. In a preliminary 
study of the directional mesh generation for viscous flow problems it has been 
found that the unstructured mesh generation scheme used here can not produce 
a good quality mesh when mesh stretching is greater than 20. For turbulent flow 
calculations highly stretched meshes, where aspect ratios are on the order of 100 


Ft I - 
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Figure 4.5 Intermediate and final mesh of NACA0012. 


to 1,000, are usually employed to resolve the shear layers. It is obvious that a 
stretching of 0(10) is not sufficient. 

In the present approach a structured triangular mesh is generated around bodies 
to_achieve high aspect ratio meshes within the boundary layer. Depending on the 
shape of the trailing edge two types of meshes are used in the present work. A 
C-type mesh is used if a blade has a wedged or cusped trailing edge and an O-type 
mesh is employed when the blade has a rounded trailing edge. 

4.2.1 O-Mesh Approach 

For blades with rounded trailing edges, once the blade boundary nodes have 
been determined the normal vector at every boundary node is defined. Radial 
running mesh lines are generated using these normal vectors. A fixed number of 
mesh points are then placed along each mesh line based on an algebraic stretching 
algorithm. That is 

ASi = Si - Si-t = V r’ -1 , i = 1, N (4.4) 

where S, indicates the normal distance from point i to the solid surface with the 
initial value So = 0, and r„ and N denote the algebraic stretching factor and total 
number of points in the normal direction, respectively. The minimum normal mesh 
scale, 6 y + , is defined as 

V = Ay wjJ i y+ (4.5) 

In the above expression A y wa u is the physical wall distance corresponding to y + = 1 
and y+ represents the allowable y + value at the first layer off the wall. 

The value of N is specified such that the aspect ratio of all cells in the outer 
layer are in the range of 10 to 20. The total thickness of the structured mesh at 
each station can be computed as: 

r N - 1 

AY^ = V ^ f 

v r 0 - 1 

The above mesh is usually referred as an O-type mesh. 


(4.6) 
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4.2.2 C-Mesh Approach 

At a wedged or cosped trailing edge the normal vectors on both upper and lower 
surfaces are different. It is appropriate to use the corresponding normal vector to 
construct mesh Unes on the upper and lower surfaces. This leads to a C-type mesh. 
The mesh point location is determined using the same procedure as for the O-mesh. 
In order to reduce the aspect ratio off the trailing edge, this mesh is extended in the 
near-body wake region. First a mesh line is constructed along the bisector vector 
of the trailing edge in the downstream direction. Then a number of grid points are 
placed on this line using the algebraic stretching distribution 


ASr. = AX 0 min{AX, 


nuuc? 


r i- l j, * = i , m 


TV - 1 


(4.7) 


where AXo is the streamwise mesh size at the trailing edge, and AX™ M 
are user-specified constants which are in the range of 4-8 and 15-20, respectively. 
The value of AX™ is specified to Unfit the growth of the streamwise mesh size at 

wake regions. 

After grid points are placed on the bisector line, mesh Unes are generated on the 
upper and lower sides of the bisector line at each node. The following distribution 
is employed to smooth the angle between the bisector vector and the upper or lower 

mesh line: 

«?=«?£ + («BI - »¥e) (4 ' 8) 

+ (*BI - *Te) jj (4 ' 9) 

In these expressions, the superscripts U and L represent the upper and lower sur- 
faces, respectively, the subscripts TE and BI denote the trailing edge and the bi- 
sector Une, respectively, the subscript i denotes the i - th mesh line, M is the total 
number of mesh Unes, and 0 is the orientation of mesh Unes. 

.Mesh points are then placed on these mesh Unes using a distribution defined as: 

A S yj = minlAFmMcVoi^-^TY}) 3 ~ ^ 


( 4 . 10 ) 
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Y 0i - Sy+ r' x \ i = 1, M (4.11) 

where Y 0 , is used to reduce the mesh aspect ratios on the wake regions. The stretch- 
ing factors, r x and r y , are in the range of 1.1-1.35. Notice that the combination of 
Eq. (4.10) and Eq. (4.11) might produce a different number of mesh points on the 
adjacent mesh lines. The use of quadrilateral meshes can only connect part of the 

mesh points. Additional triangles are therefore placed to connect the remaining 
mesh points. 

After the quadrilateral meshes are constructed a structured triangular mesh 
is generated by diagonally dividing each quadrilateral mesh cell into two triangles 
(e.g., see Figure 4.6 and Figure 4.7). The initial front list is then updated to include 
the faces of the outer layer of structured triangles in the front, while removing 
the previous fronts which connect consecutive body boundary points. With the 
new front list, the advancing front method is used to mesh the remaining domain. 
Because the algebraic stretching algorithms, Eqs. (4.4), (4.10) and (4.11), relieve 
the stretching ratio on the outer layer of the structured quadrilateral meshes, the 
advancing front method can be employed without difficulty. 



Figure 4.6 O-type structured triangular mesh. 



Figure 4.7 C-type structured triangular mesh 
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5. SOLUTION ADAPTATION ALGORITHM 

This chapter presents a solution adaptive algorithm for both inviscid and tur- 
bulent flow problems. Key elements in this algorithm are the definition of mesh 
adaptation parameters, the convergence criteria, for the flow solver and the con- 
vergence criteria for the overall adaptation loop. The mesh adaptation parameters 
are defined in terms of geometric and flow field information. For turbulent flow 
problems additional length scale information is used to define the local structured 
meshes. 

5.1 Adaptive Remeshing Algorithm 

In the present work, the flow solver and the advancing front mesh generator are 
coupled together using a solution adaptation scheme which periodically remeshes 
the solution domain to resolve the flow structure as the solution evolves. The 
basic remeshing algorithm is shown schematically in Figure 5.1. The initialization 
step generates an extremely coarse “hand-tri angulation” of the solution domain 
to be used as the initial background mesh. Mesh adaptation parameters are then 
computed at nodes on this background mesh. Since no flow field information is 
available at this point, the mesh adaptation parameters are set to constants or are 
computed based on geometric information. With this information, a coarse initial 
mesh is gen erated and t he flow solv e r is called. - • - 

The solution is marched in time on the initial mesh until a prescribed steady- 
state convergence criteria is met. Based on this initial solution, the solution adap- 
tation parameters are computed using flow field and geometric information on the 
initial mesh. The mesh is then regenerated and a new solution is computed. The 
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Figure 5.1 Flow chart of adaptive remeshing algorithm. 
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current solution is then compared to the solution on the previous mesh to check 
for overall convergence of the adaptation cycle. The adaptation cycle (outer loop 
in Figure 5.1) is repeated until the difference between the current and previous 
solutions drops below a prescribed convergence criterion. 

5.2 Mesh Adaptation Parameters 


In the advancing front method, the mesh adaptation parameters required for 
generation of a new mesh are the mesh length scale, 8, the stretching parameter, 
s, and the stretching direction denoted by or. These mesh adaptation parameters 
are computed using one or more chosen refinement parameters determined by the 
current solution and the geometry of the problem. The mesh length scale parameter, 
8, is defined in terms of two refinement parameters, one based on surface geometry 
information and the second based on flow field information. Combined, these two 
refinement parameters define the variation of 8 throughout the solution domain. 
The mesh stretching parameter, s, and the stretching direction, a, are defined 
in terms of additional flow field information such that the mesh orientation and 
stretching align with detected flow features. 


5.2.1 Geometric Refinement Parameter 


Along solid boundaries of the solution domain a refinement parameter based 
on the local surface curvature is used. This parameter has been found to be very 
important in transonic fan and compressor applications, where an accurate repre- 
sentation of the blade geometry is critical to accurately setting up the flow structure 
within the solution domain. This in turn improves the overall convergence of the 
solution adaptive method. The geometric length scale parameter, 8 a , is defined in 
terms of the local radius of curvature of the surface, R,, as 
2irRi 


(U- = 


(n.-i r 


(5.1) 


where N, is the number of points equally spaced on the arc, and (£,), is the distance 
of the arc length between any two neighboring points on the circumference of a circle 
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with radius ft. For the solutions presented in the following sections, N. is chosen 
in the range of 10 to 40 depending on the resolution desired. 

5.2.2 Flow Field Refinement Parameter 

The goal of a solution adaptive scheme is to achieve an optimal mesh on which 
the numerical error is uniformly distributed. Because the flow solution is not known 
in advance, a quantitative estimate of numerical errors must be developed. These 
errors can be estimated in many different ways. A review of various approaches is 
given by Morgan and Peraire [63]. Morgan and Peraire note that solution refine- 
ment parameters based on either first- or second-derivatives work well and are very 
economical. In the present work two sets of solution refinement parameters have 

been investigated. 

5 2.2.1 Second-Derivative Refinement Parameters 

In Reference [63] Morgan and Peraire introduce refinement parameters based 
on a local interpolation error analysis for a piecewise linear discrete approximation. 
In one dimension, the interpolation error is estimated by 

6 2 \d 2 <f>/dx 2 \ 

where 6 is the mesh spacing and <f> is some flow variable. An optimal mesh is 
determined by the requirement that the distribution of local numerical errors is 

uniform. This gives 

$ 2 1 ^ | = constant (5-2) 

dx 2 

For two-dimensional problems, the discretization errors are represented by the 
tensor of second derivatives which is given by 

d 2 <f>/dx 2 d 2 4>l dxdy 

d 2 <f>/dxdy d 2 <l>/dy 2 

This is a real symmetric tensor. For a real symmetric tensor, it can be shown that 
the eigenvalues are all real and the corresponding eigenvectors are orthogonal [4]. 
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The coordinates parallel to the eigenvectors are referred to as the principal axes. 
Morgan and Peraire invoke the orthogonality of the principal axes and introduce 
refinement parameters based on the second derivatives computed along the local 
principal directions, X\ and X 2 . This is equivalent to computing the eigenvalues of 
the local second derivative tensor. 


Ai = 


d 2 4> 


A,= 


d 2 <t> 


|Ai| > |A 2 | 


(5.3) 


dxr "‘~dxr 

To achieve an optimal two-dimensional mesh, the requirement of uniform dis 
tribution of local numerical errors is applied to all directions. This leads to 


6*|A, | = *?|A 2 | = ^ i „A mAX = constant (5.4) 

where 6 X and S 2 are the mesh spacings in the local principal directions X\ and 
X 2 , respectively, is an user-specified minimum mesh spacing, and A m , T is the 
maximum value of |A X | over the whole domain. The aspect ratio of mesh cells is 
determined as the ratio of the two mesh spacings: 

s = 6 2 f 61 = ^|Ai|/|A 2 | (5.5) 

The mesh orientation is simply defined as the angle between the principal axes and 
the Cartesian coordinates. That is 



The advantage of this approach is that the discretization error is evaluated on 
a theoretical basis and all three mesh parameters are uniquely defined. However, 
the second derivatives tend to produce large errors when there is a perturbation 
in the solution domain. For example, for transonic flow problems, a standard flow 
solver smears shock waves over several points and also produces Gibbs phenomena 
near shock wave locations. Computing second derivatives based on this solution 
gives the largest value of A j on either side of the shock wave and a small value at 
the center of the shock wave. When a new mesh is generated, small mesh cells are 




69 


constructed on either side of the shock wave and bigger mesh cells are constructed 
at the center of the shock wave. Thus, the new mesh fails to accurately resolve the 
shock wave. This problem can be alleviated by using either a high resolution flow 
solver or a refinement parameter based on lower-order derivatives. In the present 
work, the latter is used. 

5.2.2.2 First-Derivative Refinement Parameters 

For the present adaptive solution method, a refinement parameter based on the 
gradient of a certain flow property is developed. The refinement parameter, the 
absolute value of the gradient of a specified flow property, <f>, is computed at every 
node. — — 

* =i vo i= ft IF+<^ (5 ' 7) 

A quantitative local numerical error is estimated by 


& \W<t>\ = S 9i 

The variation of 6 within the solution domain is defined such that the above quantity 
is constant over the whole domain: 

6 | V <j>\ = 6 gi = 6min5m*x = constant (5.8) 


where is a user-specified minimum mesh spacing, and < 7 m*x is the maximum 
value of | V 4\ over the whole domain. The mesh orientation is computed based on 

the gradient direction. 

(5.9) 


(d<f> dA 
" = t “ Wdi) 


v< £ only is used to define two mesh parameters. Although second-order derivative 
terms can be used to estimate the mesh stretching parameter, s, this leads to the 
original problem addressed in the previous section. For complicated flow structures 
the- combined effect of the mesh stretching parameter and stretching direction can 
produce very distorted meshes. To eliminate the potential for such highly distorted 
meshes, only the mesh length scale, 6 , is used. The stretching parameter, s, is set 
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to one and the stretching direction a is taken as the direction perpendicular to the 
generating segment. Specifying s and a in this fashion results in grids composed of 
roughly equilateral triangles with varying length scales. 

_ The computed < 7 ,- is smoothed to remove numerical noise and spread the high 
gradient regions. Even with this smoothing, it is necessary to limit the allowable 
range of $r,. A cutoff level for g t is determined as follows. First a referenced value of 
the gradient of <f> is defined by the difference between the maximum and minimum 
values of <f> divided by the minimum length scale of the previous mesh. That is 


9nS = 


^min 


(5.10) 


where is the minimum length scale of the previous mesh. The cutoff value for 
g is then set at 

f 

if <7ref /jmin ^ C 2 

(5.11) 

< 7 ref> otherwise 

where C\ and C 2 are constant, typically in the range of 10-30 and 50-100, respec- 
tively. Finally, the limited value of the gradient at each point is given by 


Scut 


= < 


- _ ffmaxgcutfft f 

9* I / r (5.12) 

9m*x9c ut T (9m*x 9cu.t )9i 

This function has a maximum value of g^t and limits g in a smooth fashion. The 
mesh length scale is now determined by 


— minima*, — ^nun) (5.13) 

9i 

The minimum and maximum values of S must be specified for each mesh. The 
maximum length scale, S m defines the size of the largest triangles in the domain. 
Within the solution domain, the minimum length scale, is sequentially de- 
creased from on the initial mesh down to a value which leads to a solution 
satisfying the adaptive convergence criteria. The goal is to refine subsequent meshes 
as quickly as possible, without introducing wasted points as the solution evolves. 
Currently, this sequence is set in terms of the ratio of of the last mesh divided 
by Sjojg of the new mesh. 


[IT 1 
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5.2.2.3 Flow Feature Indicator 

Without a proper choice of the Bow feature indicator, neither the second-order 
nor the first-order error estimator can accurately and efficiently resolve the flow. 
The location of flow structures is not known in advance. An indicator related 
directly to an earlier solution must be developed. The performance of various 
indicators for inviscid transonic flows has been investigated by Dannenhoffer [21]. 
With the inclusion of viscous flow features in his results, the performance of different 
indicators for viscous flows is summarized in Table 5.1. Although a combinafon of 
different indicators have been proposed for the resolution of viscous transonic flow 
features [401, to the present work, an indicator based on the local flow speed has 

been found to be sufficient. 


5 2.2.4 Refinement Parameters for Local Structured Mesh 

The mesh spacing, 6 , computed from Eq. (5.4) or Eq. (5.13) depends on a sing 
mesh scale fw For inviscid flow calculations, the single length scale criterion is 
sufficient for resolving all flow features. Viscous flows generally involve two or more 
different length scales. For example, for turbulent transonic airfoil calculations the 
length scale required inside the boundary layer is usually on the order of 0(Re L '), 
while outside the boundary layer the length scale is on the order of the chord. 
Although one may argue that since the thickness of shock waves is on the order of 
a molecular free path, it should be reasonable to use a viscous length scale m the 
region of the shock wave. However, such a length scale increases the computing 

expense and is not required in practice. 

For the turbulent flow calculations presented here, the mesh scales for inviscid 

and viscous regions are considered separately. In the inviscid flow region the mesh 
spacing is computed using Eq. (5.13). In the boundary layer the mesh scale in the 
normal direction to the wall is based on the wall distance or wall units, y + , while the 
mesh spacing in the tangential direction to the wall is estimated using the inviscid 
mesh scale. These two mesh scales usually results in a highly stretched mesh. In 
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Table 5.1 Expected performance of different indicators for viscous transonic flow 

problems. 
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order to achieve such high aspect ratio mesh cells, a structured triangular mesh 
is employed in the present work. For the structured triangular mesh algorithm, 
Eq. (4.4), the stream wise and the normal mesh scales are specified along the wall. 
Then, a number of layers of structured mesh are generated from the wall boundary 
to the interior. The use of a structured mesh simplifies the multiple length scale 
problem in that the streamwise mesh scale specified on the wall surface is sufficient 
to define the streamwise mesh scale in the boundary layer. 

For the initial mesh the wall distance is estimated using the local skin friction 

for a turbulent flat plate flow [76]. 

Twafl = 0.0296 Re~ lf * (5-14) 

(P)V oo 

The wall distance corresponding to y + = 1 at the end of the plate is then obtained 

Ay wJ1 = 5.81ffeZ 9/W (515) 

Specifying the wall units y+ T gives the wall distance of the first structured mesh line 
off the wall: 

6„+ = Aj/waiiypr (5.16) 

where the value of 6 y + is specified in the range of 2.5-3.5 for the Baldwin-Barth 
one-equation model. The streamwise mesh scale is taken directly from the initial 
mesh scale distribution on the background mesh. From these two mesh scales, the 
structured mesh is then constructed. 

As the flow solution evolves, the local wall distance is recomputed at each bound- 
ary node using Eq. (2.14) and the minimum value of the local wall distance is used 
to estimated the new £„+ . The local streamwise mesh scale is then estimated from 
the inviscid flow information. In the inviscid flow region the Bernoulli equation 

gives 

5 + «S-° (517) 

where q is the flow speed, s is in the streamwise direction, and p is the static 
pressure. After extrapolating the inviscid flow speed from the boundary layer edge 



to the wall, the streamwise gradient of the flow speed at the wall is obtained as 


_ \LHl\ 

9 * dJ \ pqds* 


(5.18) 


where s is now replaced with the arc length scale along the wall. Employing 
Eq. (5.12) to smooth the distribution of \dqfds\ gives the streamwise mesh spacing 
along the wall. 

S H = min(<5 mai , ““^min) (5.19) 

g»i 

In the present work two additional mesh scales are used to define the surface 
mesh scale: 

S w = min{6 c , S ti , (£„+ AJ), (£„+ r^ -1 A e )} (5.20) 

where A w is the mesh aspect ratio along the wall and Ae is the mesh aspect ratio 
along the outer edge of local structured meshes. The value of A w is in the range of 
300-3,000 and A e is in the range of 10-20. From the mesh scales S y + and S w , the 
local structured mesh is regenerated. 


5.3 Convergence Criteria 

To complete the description of the present adaptive remeshing scheme the con- 
vergence criteria for the solver and adaptive remeshing cycle must be defined. Con- 
vergence for the flow solver is based on three properties of the solution: the average 
of the absolute value of the local change in pV m divided by the local time step, 
which is a measure of the convergence of the solution to steady-state; and two 
global norms defined by the difference between the highest and lowest values of the 
normal force coefficient, C/„, and the tangential force coefficient, C/ t , over the last 
50 iterations divided by reference quantities, K\ and K 2 , respectively. That is 


AC,. = [(O.U - (C7.UJ/*, 


(5.21) 


and 


AC/, — [(Cf t )niAx — (C/ t )mi n ]/f^2 


( 5 . 22 ) 
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with the reference quantities Kx and K 2 defined as 

K x = max(l, (C/Jn**), and K 2 = max(l, (C/.W) (5.23) 

The first quantity, S(pV m )/St is inversely proportional to the local time step. As 
the mesh is refined St decreases which makes this norm more restrictive as the 
mesh is adapted. The reference quantities Kx and K% are used to normalize the 
global norms. Since the flow properties are normalized with freestream values in 
the present calculations, it is appropriate to normalize the integral quantities with 
a value of 1 if the integral quantities are much less than 1. The use of windowing of 
integral quantities to define convergence of the force coefficients was developed by 
Dannenhoffer [21] for an unstructured quadrilateral based adaptive mesh scheme. 
When any one of these three coefficients drops below the user defined constant, £», 

the solution on the current mesh is converged. 

Convergence for the adaptive solution cycle is based on the change in force 
coefficients from one mesh solution to the next. When the relative change in C/ n is 
less than t 2 and the relative change in Cj, is less than £3 the adaptive mesh cycle 
is converged. The relative change is defined as 

A C f = |(C/)u* “ (C7)n*w|/max (l,max(|(C/)u. t |, |(C»m»I)) (5-24) 

The values of £i, £ 2 , and £3 are user specified constants. Specifying smaller values 
of £2 and £3 will increase the number of mesh adaptation cycles. This in turn leads 
to a smaller mesh scale and improves the accuracy of the solution. Therefore, by 
choosing these constants a user can specify the accuracy of computed solution. Note 
the value of e x should be at least one-order less than the values of c 2 and £ 3 . In the 
present work, £1 is set to 0.0001 for both inviscid and turbulent flow problems. For 
inviscid flow problems, both e 2 and £3 are set to 0.002. For turbulent flow problems, 
these values are raised to 0.005 to reduce the number of mesh adaptations. 
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6. INVISCID FLOW RESULTS 


The present adaptive solution scheme has been formulated to solve both two- 
dimensional planar flows and quasi-three-dimensional axial-, radial- and mixed-type 
flows. In the development of this approach a range of inviscid problems has been 
solved. To validate the present approach the test cases presented here include a 
two-dimensional multi-element airfoil, the Sanz supercritical compressor and sub- 
critical turbine cascades, and the Denton supersonic staggered wedge cascade. Each 
of these test cases is two dimensional and has a known analytic solution. These test 
cases cover a wide range of Mach numbers and each has different flow structures 
which must be resolved for an accurate solution. To demonstrate the adaptive solu- 
tion scheme in practical turbomachinery applications, the quasi-three-dimensional 
analysis has been used to analyze the NASA Rotor 67 low-aspect-ratio transonic 
fan and the Allison tandem blade cascade. 


6.1 Multi-Element Airfoil Case 


The first case is a model three-element airfoil operating in a high lift configura- 
tion. This case illustrates the ease with which flow problems in arbitrary multiply- 
connected regions can be solved. Reference [80] provides the airfoil geometry and 
an analytic incompressible potential flow solution for an angle of attack of 20°. In 
order to make a direct comparison with the analytic solution, this test case should 


be run with a very low freestream Mach number. Since the stability criteria for the 


present explicit scheme is inversely proportional to the speed of sound, which will 
be very large, low Mach number flows require very small time steps. Further, very 
low Mach number flows imply very small convection speeds. Both these conditions 


m i 
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result in a very slow convergence rate to the steady state solution and in turn large 
computation time. As a compromise between the computation time and the com- 
pressibility effects, this test case was run with freestream Mach number of 0.125. 
The flow conditions and airfoil characteristics are summarized in Table 6.1. The 
far field and near field meshes are shown in Figures 6.1 and 6.2, respectively. The 
far field boundary is set at a radius of 5 chords from the airfoil. At this boundary 

the vortex far field boundary condition [83] is used. 

The initial mesh is refined near the airfoil through the specification of the length 
scale parameter along the surface (N, = 40). In the far field the mesh is uniform 
with a length scale of 0.5 chord. Within the solution domain the minimum 

length scale is sequentially decreased, from on the initial mesh, down to 
a value which leads to a solution satisfying the adaptive convergence criteria. The 
refinement sequence is set in terms of the ratio of of the last mesh divided by 
of the new mesh. For the airfoil case presented here this ratio equals 32 for 
the first adaptation, 4 for the second and the third adaptations and 2 for the last 
refinement. The sequence of solution adapted meshes are shown in Figures 6.2 to 
6.6. Mesh statistics and force coefficient information for each mesh are summarized 

in Table 6.2. 

Mach number contours for the final solution are shown in Figure 6.7. Note 
that the maximum surface Mach number on the leading edge slat is about 0.6, 
which is well out of the range of what can be considered incompressible. Figure 6.8 
compares the converged surface pressure coefficient distribution over each element 


Table 6.1 Flow conditions and airfoil characteristics for the three-element airfoil. 


Moo 

a 

£«ut 


C^exact 

0.125 

20.0° 

45.0° 

15.0° 

5.136 

Note: 

Incompressible analytic solution 


is available in Reference [80]. 
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with the incompressible analytic solution. Although the flow structure of this test 
case involves high gradients in flow properties, they are only concentrated in small 
regions where the surface curvature is high. Therefore, plots of intermediate mesh 
solutions are almost exactly the same. Most of the error in the initial solution is 
concentrated in the region of the suction peak on the slat. 

To illustrate the importance of compressibility in this case, a compressibility 
correction to the analytic solution using the Karman-Tsien rule is considered in 
Figure 6.9. The compressible presure coefficient using the Karman-Tsien rule is 
given as: 



where C pfi is the incompressible pressure coefficient, is the freestream Mach 
number, and C p is the corrected pressure coefficient. Since the high Mach number 
flows are concentrated in the suction surface of the slat portion and the leading 
edge suction surface of the airfoil, the corrected surface pressure coefficients are 
lower than the incompressible values in these regions. Otherwise, the corrected 
surface pressure coefficients are essentially the same as the incompressible values. 
Figure 6.9 shows a slight difference between the computed pressure coefficient and 
the corrected pressure coefficient in the high Mach number flow regions. This is 
due to the extremely high flow gradients in these regions. These regions can be 
accurately resolved with further mesh adaptation. 

Figures 6.10, 6.11, and 6.12 show the surface total pressure errors for the initial, 
3 adapted, and final meshes, respectively. The total pressure error is defined as 
the difference between local total pressure and the upstream value. That is, 

p _ 1 _ gv 

terror p 

*Tq 

Since this flow is isentropic the total pressure error should be zero. These plots 
clearly show a reduction in error as the mesh is refined in the slat region. There is 
a 40 % reduction in the peak total pressure error in the final mesh solution relative 
to the 3rd adapted mesh solution. The final mesh length scale in the slat leading 
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edge region is about 60 % of the 3 rd adapted mesh length scale in the same region. 
This indicates that the present scheme is first-order. 

The computed lift coefficient is 3 % less than the Karman-Tsien corrected lift 
coefficient which is about 5.279. It is believed that part of this error is due to the 
location of the far field boundary. Based on the results of Usab [83] and Dannenhof- 
fer [21], the lift coefficient increases about 1.5 to 2.0 % as the far field boundary is 
moved from 5 chord to 50 chord. It should also be noted that for low Mach number 
flows, the total pressure error has a significant effect on the pressure distribution. 
Even though the total pressure error in the final mesh is very small, if all the total 
pressure error is assumed to reduce the surface static pressure this would lead to a 
5 % error in the predicted lift coefficient. Although a direct relation between the 
lift coefficient and total pressure error is not known, it is reasonable to expect a 1 
to 2 % reduction in the lift coefficient due to this total pressure error. 

Finally, the convergence histories for the lift coefficient and average residual, 
\S(pu)/6t\, are shown in Figures 6.13 and 6.14, respectively. Note the reduction in 
the convergence rate of the flow solver on the latter meshes. This is due to the 
decreasing mesh length scale in the fine mesh region where the explicit time step is 

very small. 

Table 6.2 Three-element airfoil: Mao = 0.125, a = 20°, and (C'Oincomp = 5.136. 


MESH 

0 

1 

2 

3 

4 

nodes 

2,034 

2,092 

2,911 

4,643 

8,506 

elements 

3,769 

3,888 

5,517 

8,927 

16,567 

Ci 

4.9892 

4.9847 

5.0630 

5.1061 

5.1165 

Ci 

0.0570 

0.0444 

0.0286 

0.0201 

0.0145 

Sqm/ chord 

0.5 

0.5 

0.5 

0.5 

0.5 

^miT j f Sfgun 

1 

32 

128 

512 

1,024 


Total CPU: 1,307 sec in Cray-YMP. 
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Figure 6.2 Multi-element airfoil initial mesh: 2,034 nodes and 3,769 elements. 
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Figure 6.3 Multi-element airfoil 1“ adapted mesh: 2,092 nodes and 3,888 

elements. 



Figure 6.4 Multi-element airfoil 2 nd adapted mesh: 2,911 nodes and 5,517 

elements. 


Figure 6.5 Multi-element airfoil 3^ adapted mesh: 4,643 nodes and 8,927 

elements. 



Figure 6.6 Multi-element airfoil 4 th (final) adapted mesh: 8,506 nodes and 16,567 

elements. 


Figure 6.7 Multi-element airfoil Mach number contours on the final mesh 
cmin = 0.0, cmax = 0.60, and inc = 0.02. 


Figure 6.8 Multi-element airfoil surface pressure coefficient for the final mesh 
solid line - numerical solution and symbol - analytic solution. 
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6.2 Sanz Supercritical Compressor Cascade 

This case is a shock free supercritical compressor cascade designed by Sanz [75] 
using the hodograph method. This case is described in the AGARD report [25] 
as test case A/CA-1. The hodograph method used to design this cascade results 
in a blade section with a nonclosed trailing edge. As noted in Reference [25], any 
attempt to close the blade trailing edge changes the blade section, resulting in 
a different solution. The flow conditions and blade characteristics are listed in 
Table 6.3. Afexit is the specified isentropic exit Mach number which sets the ratio 
of the exit static pressure to the inlet total pressure. This ratio is given as 

^ 7=( 1 + — ( 6 - 2 ) 

The ratio is 0.8177 for the current value of A/ exit . The inflow and outflow boundaries 
for the mesh are placed one half an axial chord from the blade leading and trailing 
respectively. At the inflow boundary total pressure, total temperature, and 
absolute flow angle are specified. At the outflow boundary the exit pressure is set. 
Flow tangency is enforced on the blade surface and periodicity is imposed at the 
upper and lower boundaries of the domain. At the open trailing edge a uniform jet 
flow condition is specified. 

The sequence of solution adapted meshes is shown in Figures 6.15 to 6.16. 
Because the flow structure for this test case does not involve high gradients it only 

Table 6.3 Flow conditions and blade characteristics for the Sanz supercritical 

compressor cascade. 



•Wexit 

An 

Axit 

Pitch/chord 

0.711 

0.544 

30.81° 

-0.35° 

1.034 


Note: • Hodograph solution is available 


in Reference [25]. 

• Nonclosed trailing edge present 
in the physical profile. 
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Table 6.4 Sana supercritical compressor cascade: M m - 0.711 and - 30.81 


MESH 

0 

1 

nodes 

820 

947 

elements 

1,469 

1,717 

Cfn 

0.9044 

0.9039 

Cf, 

-0.2817 

-0.2825 

$m»x/ pitch 

0.095 

0.095 

&mMxj ^min 

1 

4 

Total CPU: 34 sec in Cray-YMP. 


takes one remesh cycle to converge. The mesh and loading coefficient information 
are given in Table 6.4. Mach number contours for the final solution are shown 
in Figure 6.17. A comparison of the converged Mach number distribution over the 
blade with the hodograph solution is shown in Figure 6.18. The agreement between 
the numerical and analytic solutions is excellent except in a small region near the 
sonic point on the suction surface. Note that the supercritical blade geometry is an 
isolated shock free design and that very small variations in the geometry will lead 
to the formation of shock waves. Since the blade coordinates given in Reference [25] 
have gaps in the sonic regions of the blade, the blade geometry in these regions is 
defined through interpolation. The difference between the computed solution and 
the analytic solution is the result of the difference between the interpolated and 
the unknown correct blade geometry in these regions. Convergence histories for the 
normal force coefficient and average residual are shown in Figures 6.19 and 6.20, 

respectively. 
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Figure 6.15 Sanz supercritical compressor cascade initial mesh: 820 nodes and 

1,469 elements. 



Figure 6.16 Sanz supercritical compressor cascade 1 st (final) adapted mesh: 947 

nodes and 1,717 elements. 



Figure 6.17 Sanz supercritical compressor cascade Mach number contours on the 
final mesh: cmin = 0.0, cmax = 1.30, and inc = 0.05. 



Figure 6.18 Sanz supercritical compressor cascade surface Mach number for the 
final mesh: solid line - numerical solution and symbol - analytic solution. 
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6.3 Sanz Subcritical Turbine Cascade 

This subcritical turbine cascade, AGARD test case A/CA-3 [25], was also de- 
signed using a hodograph method. The blade section has a nonclosed trailing edge. 
The flow conditions and blade characteristics are tabulated in Table 6.5. The up- 
stream and downstream mesh boundaries are placed one half an axial chord from 
the blade edges. The numerical boundary conditions are the same as those used in 
the supercritical compressor cascade case. The ratio of the exit static pressure to 
the inlet total pressure is 0.6788 for the present value of M exit . 

Figures 6.21 to 6.24 present the sequence of solution adapted meshes. Since the 
gradients are higher in this case due to the rapid acceleration of the flow through 
the turbine blade passage, it takes three mesh adaptation cycles to converge. Mesh 
properties and force coefficients for this case are summarized in Table 6.6. Mach 
number contours for the final solution are shown in Figure 6.25. Figures 6.26 and 
6.27 compare the computed surface Mach number distribution with the hodograph 
solution for the initial and final meshes, respectively. These plots clearly show the 
improvement in the numerical solution as the mesh is refined. 

For inviscid flows, the total pressure error is also a good indicator of numerical 
error in the solution scheme. Since this flow is isentropic the total pressure error 
should be zero. The plots of surface total pressure error shown in Figures 6.28 

Table 6.5 Flow conditions and blade characteristics for the Sanz subcritical 

turbine cascade. 



Merit 

Am 

mm 

Pitch/chord 

0.343 

0.765 

36 . 0 ° 

-57.35° 

0.566 

Note: 

• Hodograph solution is available 

in Reference [25]. 

• Nonclosed trailing edge present 

in the physical profile. 


IT I 
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and 6.29 demonstrate the improved accuracy of numerical solution as the mesh is 
refined. The maximum surface total pressure error is decreased from 4 % on the 
initial mesh down to 2 % on the final mesh. Convergence histories for the normal 
force coefficient and average residual are given in Figures 6.30 and 6.31, respectively. 


Table 6.6 Sanz subcritical turbine cascade: M m = 0.343 and An - 36.0 



Figure 6.21 Sanz subcritical turbine cascade initial mesh: 722 nodes and 1,244 

elements. 
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Figure 6.22 Sanz subcritical turbine cascade 1 st adapted mesh: 1,018 nodes and 

1,799 elements. 



Figure 6.23 Sanz subcritical turbine cascade 2 nd adapted mesh: 1,620 nodes and 

2,962 elements. 
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Figure 6.24 Sanz subcritical turbine cascade 3 rd (final) adapted mesh: 2,336 

nodes and 4,367 elements. 



Figure 6.25 Sanz subcritical turbine cascade Mach number contours on the final 
mesh: cmin = 0.0, cmax = 0.80, and inc = 0.05. 
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6.4 Denton Supersonic Staggered Wedge Cascade 

The Denton wedge cascade, AGARD test case A/CA-2 [25], is a compressor 
cascade operating in a fully supersonic flow. This staggered wedge cascade has 
been carefully designed so that the reflected obUque shock wave is exactly cancelled 
resulting in an uniform flow between the two parallel surfaces. In order to close 
the cascade profile, another wedge is introduced on the pressure side resulting in 
an expansion off the downstream comer. Figure 6.38 presents the exact solution 
derived from shock/expansion theory. This test case illustrates the capability of the 
present solution adaptive scheme in capturing complicated shock wave structures. 
The flow conditions and blade characteristics are listed in Table 6.7. 

The initial mesh is shown in Figure 6.32. The upstream and downstream mesh 
boundaries are located one half an axial chord from the blade leading and trailing 
edges, respectively. At the inflow boundary total pressure, total temperature, and 
absolute flow angle are specified as boundary conditions. At the outflow boundary 
the nonuniform static pressure boundary condition is imposed (see section 3.5.4). 

Figures 6.32 to 6.35 present the mesh sequence during the solution adaptation 
process. The mesh characteristics are summarized in Table 6.8. Comparing the 
initial mesh Mach number contours shown in Figure 6.36 with the final mesh Mach 
number contours in Figure 6.37 shows an improved shock wave resolution for the 

Table 6.7 Flow conditions and blade characteristics for the Denton supersonic 

staggered wedge cascade. 


Min 

Mexit 

An 

/^exit 

Pitch/chord 

1.6 

1.401 

60.0° 

60.0° 

0.47985 


Note: • Analytical solution is available 


in Reference [25]. 

• Nonuniform pressure present 
at the exit. 
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Table 6.8 Denton supersonic staggered wedge cascade: - 160, An - 60.0° 

(C /t W = 0.03774, and (<?/„)««* - -0.0654. 


MESH 


i 

2 

3 

nodes 

1,231 

1,982 

3,625 

7,526 

elements 

2,196 

3,562 

6,762 

14,440 

c Jm 

0.0411 

0.0375 

0.0384 

0.0382 

Cj, 

-0.0722 

-0.0647 

-0.0667 

-0.0664 

SgtMxl pitch 

urn 

0.08 

0.08 

0.08 

^upar/ ^min 

i 

3.1 

12.5 

25 


Total CPU: 880 sec in Cray-YMP- 


adapted solution. Mach number contours for the final and theoretical solutions are 
shown in Figures 6.37 and 6.38, respectively. Even though the shock wave structure 
is complex, the remeshing procedure accurately resolves the flow. The computed 
surface Mach number distributions on the initial and final meshes are compared 
with the theoretical solution in Figures 6.39 and 6.40, respectively. These plots 
show a great improvement in the accuracy of numerical results. Since the present 
scheme is not a monotonic scheme, overshoots are shown in the shock wave regions. 
Comparing the computed loading coefficients C/„ and Cj t to the analytical values 
in Table 6.8 shows a reduction in the error as the mesh is refined with the exception 
of the l rt adapted mesh solution. This indicates that it is difficult to determine the 
improvement of solution based on the loading coefficients only. Since the loading 
coefficients are integral quantities, numerical errors may cancel after integrating 

over the blade surface. 

Convergence histories are given in Figures 6.41 and 6.42. It is noted that the 
CPU time this case requires is five times longer than that of the Sanz subcntical 
turbine cascade case. This is due to the increasing number of mesh points and 
decreasing mesh length scale as the mesh is adapted to the solution. 
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Figure 6.32 


Denton supersonic staggered wedge cascade initial mesh: 1,231 nodes 
and 2,196 elements. 


inn it 
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Figure 6.33 Denton supersonic staggered wedge cascade 1“ adapted mesh: 1,982 

nodes and 3,562 elements. 
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Figure 6.34 Denton supersonic staggered wedge cascade 2 
® nodes and 6,762 elements. 


iirir 


adapted mesh: 3,625 
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Figure 6.35 Denton supersonic staggered wedge cascade 3 rd (final) adapted mesh: 

7,562 nodes and 14,440 elements. 
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ade Mach number contours 
52, and inc = 0.03. 
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Figure 6.37 Denton supersonic staggered wedge cascade Mach number contours 
on the final mesh: cmin = 1.17, cmax = 1.62, and inc = 0.03. 
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analytic solution (Denton tt al. [25]). 




109 



Figure 6.39 Denton supersonic staggered wedge cascade surface Mach number for 
the initial mesh: solid line - numerical solution and symbol - analytic solution. 
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Figure 6.40 Denton supersonic staggered wedge cascade surface Mach number for 
the final mesh: solid line - numerical solution and symbol - analytic solution. 
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Figure 6.41 Denton supersonic staggered wedge cascade: C/ n verses iteration. 



Figure 6.42 Denton supersonic staggered wedge cascade: average \8(pu)jSt\ verses 

iteration. 
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6.5 NASA Rotor 67 Transonic Fan 

NASA rotor 67 is a low- aspect-ratio transonic axial-flow fan rotor. Laser anemome- 
ter surveys of the flowfield were made for operating conditions near peak efficiency 
and near stall by Strazisar et d. [79]. The fan tip relative Mach number is 1.38. The 
experiments were conducted in a rotor-only configuration making them good test 
cases for the three-dimensional flow solver. Only the near peak efficiency test con- 
dition will be considered here. Since the flow in this machine is three-dimensional, 
axisymmetric through- flow information is required as input to the present quasi- 
three-dimensional analysis. The streamsurface location and thickness data used in 
the present calculations were obtained from Reference [60]. It is important to note 
that the accuracy of the present results depends on the accuracy of the streamsur- 
face data, which is difficult to validate. This test case is important in illustrating 
how the adaptive solution procedure resolves realistic flow structures. 

The adaptive mesh solution is computed for the 30 %-span streamsurface (mea- 
sured from the shroud). In this case the inlet boundary for the mesh is placed 
three quarters an axial chord from the blade leading edge, and the exit boundary 
is located half an axial chord from the blade trailing edge. At the inflow boundary 
total pressure, total temperature, and absolute flow angle are specified as boundary 
conditions. At the outflow boundary the exit pressure is imposed. Since there is 
high turning involved in this flow, it is not possible to start a solution with constant 
initial conditions. Therefore, the quasi-one-dimensional solution is used to provide 
a smooth variation of initial flow conditions (see section 3.4). The inflow and out- 
flow condition as determined from the through-flow analysis are summarized in 

Table 6.9. 

This span station is particularly interesting because the upstream relative Mach 
number is 1.20. With a blunt leading edge, a bow shock wave stands away from the 
blade leading edge. To accurately predict this flow the leading edge must be well 
resolved. Referring to Figures 6.50 and 6.52 it is clear that a refinement parameter 
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Table 6.9 NASA Rotor 67 far field flow conditions. 


Case 

P* i.(P s} ) 

n. m 

rv$ 

Pexit/Ftu 

30 % span 

14.7 

Cn 

h- * 

po 

o 

0 

1.21 

50 % span 

14.7 

518.7° 

0 

1.25 

70 % span 

14.7 

518.7° 

0 

1.16 


based on surface curvature does a good job in resolving the leading edge region. 
Even on the initial mesh the bow shock wave is located in the proper location. 
Without this refinement the adaptive procedure improperly locates the bow shock 
wave. This changes the flow conditions within the blade passage, resulting in a very 
poor initial solution. It then takes many mesh adaptation cycles to converge to the 
correct solution. The present adaptation criteria leads to a converged solution after 
three remeshes (see Figures 6.43 to 6.46). Mesh properties and force coefficients 
are summarized in Table 6.10. 

The relative Mach number contours for the final mesh are shown in Figure 6.47. 
The bow shock reflection on the blade suction surface is not well resolved, because 
the shock strength is much weaker than the normal shock. Figure 6.48 presents 
contours of the experimentally measured relative Mach number for the 30 %-span 
station. Even though a very large amount of experimental data was taken, the 
experimental data used in making these contour plots are still very sparse. This fact 
makes it very difficult to evaluate the overall accuracy of the computed solution. 
The compute results show a good agreement with the experimental data in the 
upstream of the blade leading edge. Note the experimental data shows a passage 
shock wave located further upstream than the shock wave location computed here. 
Boundary layer blockage plays an important role in determining the location of the 
passage shock wave. Since this is an inviscid solution the shock wave is located 
downstream of the correct experimental location. 

Figures 6.51 and 6.52 show a blowup of the leading edge region of the blade. 
There are about 20 mesh points along the half circle of the leading edge. The bow 


HTTP 
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shock wave is well resolved. Plots of the surface pressure coefficient are shown in 
Figures 6.53 and 6.54 for the initial and final meshes, respectively. Convergence 
histories are shown in Figures 6.55 and 6.56. The normal force coefficient (force in 
the ^-direction) converges rapidly to a constant value. The adaptation cycle has 
been stopped at the 3 rd adapted mesh, since further refinement will lead to a large 
number of mesh points which in turn increases the computation time. The relative 
change in the force coefficient between last two meshes is within 1 %. 

Table 6.10 NASA Rotor 67 operating at peak efficiency: 30 % span station. 


MESH 

0 

1 

2 

3 

nodes 

1,327 

2,031 

5,138 

7,245 

elements 

1 

2,375 

3,700 

9,877 

14,074 

c fn 11 

1.2149 

1.1919 

1.1950 

1.2073 

\c h 

-1.6292 

-1.5763 

-1.5931 

-1.6270 

&maxl pitch 

0.067 

0.067 

0.067 

0.067 

^mxY j ^min 

1 

4 

16 

32 

Total CPU: 

553 sec in Cray-YMP. 
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Figure 6.43 NASA Rotor 67 30 % span from shroud initial mesh: 1,327 nodes 

and 2,375 elements. 
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Figure 6.44 NASA Rotor 67 30 % span from shroud l* 1 adapted mesh: 2,031 

nodes and 3,700 elements. 
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Figure 6.45 NASA Rotor 67 30 % span from shroud 2 nd adapted mesh: 5,138 

nodes and 9,877 elements. 
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Figure 6.46 NASA Rotor 67 30 % span from shroud 3 rd (final) adapted mesh: 

7,245 nodes and 14,074 elements. 


Figure 6 47 NASA Rotor 67 30 % span from shroud relative Mach number 
contours on the final mesh: cmin = 0.0, cmax = 1.70, and inc = 0.05. 
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Figure 6.48 NASA Rotor 67 30 % span from shroud relative Mach number 
contours: experimental data (Strazisar, et al. [79]). 
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Figure 6.49 NASA Rotor 67 30 % span from shroud: leading edge blowup of the 

initial mesh. 



' Finite 6.50 NASA Rotor 67 30 % span from shroud relative Mach number 
contours (cmin = 0.0, cmax = 1.70, and inc = 0.05): initial mesh leading edge 

region. 



Figure 6.51 NASA Rotor 67 30 % span from shroud: leading edge blowup of the 

final mesh. 



Figure 6.52 NASA Rotor 67 30 % span from shroud relative Mach number 
contours (cmin = 0.0, cmax = 1.70, and inc = 0.05): final mesh leading edge 

region. 
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Figure 6.55 NASA Rotor 67 30 % span from shroud: C Jn verses iteration. 



Figure 6.56 NASA Rotor 67 30 % span from shroud: average \S(pV m )/8t\ verses 

iteration. 


II FT 



123 


6.6 Allison Tandem Blade Cascade 

The Allison tandem blade cascade is an axial-flow compressor stator configura- 
tion. This is a preliminary design and there are no experimental data available for 
this case. This case demonstrates the ease with which arbitrary multiply-connected 
regions ran be solved. Quasi-three-dimensional through-flow information is used in 
the present quasi-three-dimensional analysis. The streamsurface location and thick- 
ness data used in the present calculation axe obtained from Reference [60]. 

The adaptive mesh solution is computed for the 70 %-span streamsurface. The 
far field and near field meshes are shown in Figures 6.57 and 6.58, respectively. In 
this case the far field boundary for the mesh is placed one axial chord from the 
blade edges. Numerical boundary conditions are the same as those used for the 
NASA Rotor 67 case. Uniform initial flow conditions are used. The inflow and 
outflow conditions determined from the through-flow analysis are summarized in 
Table 6.11. 

Table 6.11 Allison tandem blade cascade far field flow conditions. 


Case 

P*u( P si ) 

T,„ (R) 

rv$ 

Pexit/A. 

30 % span 

74.5 

872.7° 

0 

0.791 

50 % span 

74.5 

872.7° 

0 

0.793 

70 % span 

74.5 

872.7° 

0 

0.796 


The remesh sequence is shown in Figures 6.57 to 6.60. These plots illustrate that 
flow problems of arbitrary multiply-connected regions are no more difficult to solve 
than flow over a single blade. Mesh statistics and loading coefficient information 
for this case are given in Table 6.12. This case takes two adaptation cycles to 
converge. Mach number contours for the final solution are shown in Figure 6.61. 
The inviscid solution shows shock waves in both leading edge regions. This is 
due to the rapid acceleration of the flow around a small finite radius leading edge. 
Viscous effects will slow the acceleration and may even result in a laminar leading 
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edge separation bubble. These shock waves will not- be present in a viscous flow 
calculation. Figures 6.62 to 6.65 show a blowup of the leading edge region and the 
gap region between the tandem blades. These plots show a smooth variation of 
mesh length scale in the gap region. This kind of mesh resolution is very difficult, if 
not impossible, to achieve using a structured mesh flow solver. Figures 6.66 and 6.67 
show velocity vectors for the leading edge and the gap regions, respectively. Due 
to the high flow incidence relative to the second blade there is a rapid acceleration 
of the flow through the gap region. 

The plots of the surface pressure coefficient are shown in Figures 6.68 and 6.69 
for the initial and final meshes, respectively. The spikes in the leading edge regions 
are caused by the rapid acceleration of the flow. Convergence histories are shown 
in Figures 6.70 and 6.71. The complete procedure converged in 3,000 iterations. 

Table 6.12 Allison tandem blade cascade: 70 % span station. 


MESH 

0 

1 

2 

nodes 

1,822 

2,415 

4,361 

elements 

3,292 

4,445 

8,279 

<?/« 

0.8141 

0.7802 

0.7885 

Cf, 

-0.3674 

-0.3689 

-0.3699 

&max/ pitch. 

0.096 

0.096 

0.096 

Smxx 1 &m\n 

1 

4 

16 

Total CPU: 478 sec in Cray-YMP. 
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Figure 6.57 


Allison tandem blade cascade 70 % span from shroud initial mesh: 
far field view. 


Figure 6.58 Allison tandem blade cascade 70 % span from shroud initial mesh: 

1,822 nodes and 3,292 elements. 



Figure 6.59 Allison tandem blade cascade 70 % span from shroud 1 
mesh: 2,415 nodes and 4,445 elements. 


st adapted 
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Figure 6.60 Allison tandem blade cascade 70 % span from shroud 2 nd (final) 
adapted mesh: 4,361 nodes and 8,279 elements. 



Figure 6.61 Allison tandem blade cascade 70 % span from shroud Mach number 
contours on the final mesh: cmin = 0.0, cmax = 1.60, and inc = 0.05. 







Figure 6.64 Allison tandem blade cascade 70 % span from shroud: tandem blade 

gap blowup of the final mesh. 



Figure 6.65 Allison tandem blade cascade 70 % span from shroud Mach number 
contours (cmin = 0.0, cmax = 1.60, and inc = 0.05): final mesh tandem blade gap 

region. 
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Figure 6.66 Allison tandem blade cascade 70 % span from shroud velocity vector: 

final mesh leading edge region. 



Figure 6.67 Allison tandem blade cascade 70 % span from shroud velocity vector: 

final mesh tandem blade gap region. 
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ITERATION 


Figure 6.70 Allison tandem blade cascade 70 % span from shroud: C/„ verses 

iteration. 



Figure 6.71 Allison tandem blade cascade 70 % span from shroud: average 

\8(pV m )/St\ verses iteration. 




133 


7. TURBULENT FLOW RESULTS 

In this chapter the solution adaptive scheme is applied to turbulent flow prob- 
lems. The following turbulent flow problems are considered: a subsonic flat plate; 
a RAE2822 airfoil; a NLR two-element airfoil; the VKI LS82-1 turbine cascade, 
and the Allison tandem blade cascade. Since experimental data for the first four 
cases is available, these test cases allow a validation of the present solution adaptive 
scheme. Numerical results for the Allison tandem blade cascade are presented to 
demonstrate the advantage of an adaptive unstructured turbulent flow solver for 
new advanced turbomachinery blade designs. The objective here is to illustrate 
the capability of the present solution adaptive scheme in resolving complex flow 
structures. Thus, the convergence criteria are set at a higher level to reduce the 
computing work. 

7.1 Subsonic Flat Plate 

The subsonic flat plate flow case is a standard test case for the validation of 
turbulent flow calculations. Because all the turbulence models were calibrated us- 
ing the flat plate solution, an accurate prediction of this flow problem indicates a 
correct implementation of the turbulence model equations. The purpose of this test 
case here is to validate the numerical implementation of the present turbulent flow 
solver, therefore, adaptive remeshing is not performed. The triangular mesh used 
to compute the subsonic flat plate turbulent boundary layer solution is generated 
from a 72x80 structured quadrilateral mesh. The mesh, shown in Figure 7.1, is 
algebraically packed near the leading edge and solid surface with minimum incre- 
ments of Aimm = 0.0003 and At/ ,,,;,, = 0.000003. The maximum aspect ratio in 
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this mesh is about 40,000. The inlet boundary is placed at one half of the plate 
length from the leading edge and the exit boundary is located at the end of flat 
plate. Numerical boundary conditions for the inlet and exit boundaries are the 
same as those used in the inviscid cascade cases. The freestream static pressure is 
prescribed along the exit plane. A no-slip condition with adiabatic wall and zero 
normal pressure gradient are imposed along the solid surface. Flow tangency is en- 
forced along the symmetric boundary which lies along the lower boundary upstream 
of the plate. The computed results are compared to the experimental data given in 
the 1968 AFOSR-IFP-Stanford Conference [89] for the Bladwin-Barth one-equation 
and the Chien k - c turbulence models. The flow conditions for this case are given 

in Table 7*1. 

In structured mesh approaches, it has been found that the accuracy of the 
viscous solutions is very sensitive to the effect of artificial dissipation. In order to 
determine a proper scaling of the artificial dissipation for the present unstructured 
flow solver, three different scaling of the artificial dissipation terms are investigated. 
The first scaling is the standard modified eigenvalue scaling, Eq. (3.22). The second 
scaling combines a local velocity scaling with the modified eigenvalue scaling (see 
Eq. (3.24)). Note that a local velocity scaling based on q* was also investigated in 
the present work. However, this scaling does not produce a stable solution. In the 
last scaling, the artificial dissipation is switched off if the local velocity is less than 
a specified cut-off flow speed, U mt . This scaling is a modified velocity scaling which 

is described as / / \ \ a \ 

*-« 4 + (£)) (71) 


Table 7.1 Flow conditions for the subsonic flat plate. 


Moo 

a 

Poo (psi) 

Too (° K) 

Rcl 

LrrS (m) 

0.25 

o 

O 

o 

14.7 

293.0 

2.2xl0 7 

10 

Note: 

Experimental data are available in Reference [89]. 
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where 

Q = min {l,max {o, (? 2 - C^»t)/(^« “ ^t)}} 

In the present calculation Ucat is taken as 50 % of freestream velocity. This scaling 
does not guarantee a stable solution. In the present work, the scaling, Eq. (7.1), is 
only used for the assessment of the accuracy of the turbulent flat plate flow solution. 
In order to achieve a stable solution, this scaling is vised after the flat plate solution 
is obtained with the local velocity scaling, Eq. (3.24). 

Figures 7.2 and 7.3 present the log-law velocity profile predicted by the Baldwin- 
Barth one-equation and the Chien low Reynolds number k - e turbulence models 
at Re x = 1.13xl0 7 . These plots clearly show that the eigenvalue scaling does 
not produce a satisfactory log-law velocity profile. With local velocity scaling the 
results are great im proved. Specifying zero artificial dissipation in the viscous layer 
further improves the predicted results for both turbulence models. Figures 7.4 and 
7.5 show the nondimensional velocity profile predictions for each model at x=4.987 
m. The nondimensional velocity is defined as u/U 0 o, the ratio of local velocity to 
the freestr eam velocity. Results are compared to experimental data. Results for 
the velocity scaling and the switch-off approach agree quite well with experimental 
data. Figures 7.6 and 7.7 present the coefficient of friction prediction along the 
plate for both turbulence models. Again the local velocity scaling results are better 
than the modified eigenvalue scaling. Without adding artificial dissipation in the 
viscous layer the friction coefficient distribution is accurately predicted. 

In the present turbulent flat plate flow case specifying zero artificial dissipa- 
tion in the viscous layer produces an accurate prediction. However, for general flow 
problems such an approach can lead to numerical stability problems. In the remain- 
ing cases the local velocity scaling artificial dissipation is used. Both Baldwin-Barth 
one-equation and Chien low Reynolds number k — e turbulence models have been 
successfully demonstrated on the turbulent flat plate flow. It has been recognized 
that Chien’s low Reynolds number k - e model requires a finer wall spacing (y + « 1) 
than the Baldwin-Barth one-equation model (y + « 3.5). In addition, for the low 
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Reynolds number k — c model it is very difficult to obtain a stable solution without 
a proper initialization of the turbulence quantities. In structured mesh schemes 
k and c are usually initialized using an algebraic turbulence model solution [26]. 
However, such an initialization is not practical in an unstructured mesh scheme. 
Since the Bladwin-Barth one-equation is more robust and requires less mesh reso- 
lution in the near-wall regions, only the Baldwin-Barth one-equation model is used 
in the remaining turbulent flow problems. 


T 
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Figure 7.1 Triangular mesh for flat plate flow case: 5,760 nodes and 11,218 
elements (10:1 magnification in y-direction). 



Figure 7.2 Flat plate log-law velocity profile at Re x — 1.13xl0 7 obtained with 
Baldwin-Barth one-equation turbulence model for different artificial dissipations: 
□ - eigenvalue scaling, © - local velocity scaling, A - zero artificial dissipaton in 
the bo un dary layer and solid line - standard log-law profile 
u+ = 5.85 log 10 y + + 5.56. 
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Figure 7.3 Flat plate log-law velocity profile at Re x = 1.13xl0 7 obtained with 
Chien’s low Reynolds number k — e turbulence model for different artificial 
dissipations: □ - eigenvalue scaling, © - local velocity scaling, A - zero artificial 
dissipaton in the boundary layer and solid line - standard log-law profile 

u + = 5.85 log 10 y + + 5.56. 
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Figure 7.4 Flat plate nondimensional velocity profile at x=4.987 m obtained with 
Baldwin-Barth one-equation turbulence model for different artificial dissipations: 
solid line - eigenvalue scaling, dashed line - local velocity scaling, short dashed line 
- zero artificial dissipaton in the boundary layer and symbol - experimental data. 
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Figure 7.5 Flat plate nondimensional velocity profile at x=4.987 m obtained with 
Chien’s low Reynolds number k — e turbulence model for different artificial 
dissipations: solid line - eigenvalue scaling, dashed line - local velocity scaling, 
short dashed line - zero artificial dissipaton in the boundary layer and symbol - 

experimental data. 
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x 

Figure 7.6 Flat plate surface friction coefficient obtained with Baldwin-Barth 
one-equation turbulence model for different artificial dissipations: solid line - 
eigenvalue scaling, dashed line - local velocity scaling, short dashed line - zero 
artificial dissipaton in the boundary layer and symbol - experimental data. 
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Figure 7.7 Flat plate surface friction coefficient obtained with Chien’s low 
Reynolds number k — c turbulence model for different artificial dissipations: solid 
line - eigenvalue scaling, dashed line - local velocity scaling, short dashed line - 
zero artificial dissipaton in the boundary layer and symbol - experimental data. 
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7.2 RAE2822 Airfoil 

The first test case presented for the solution adaptive scheme is the 1987 VTA 
Workshop test case Bl. This case is a RAE2822 airfoil operating at Moo = 0.725, 
a = 2.92°, and Rei = 6.5xl0 6 . This flow involves a strong shock wave at 55 % 
chord on the upper surface. The lift coefficient in this case depends strongly on 
the predicted shock wave location. This makes a good resolution of the shock wave 
very important. Experimental data is given by Cook, McDonald, and Firmin [20]. 
Note that a transition trip was used near the leading edge suction surface in the 
experiment. The flow conditions used for the present calculation are listed in Ta- 
ble 7.2. For this test case, each participant of the 1987 VTA Workshop used a 
different corrected angle of attack in his numerical calculations to obtain a better 
agreement with experimental data. The range of corrected angles varied from 2.3 to 
2.92 in the workshop [32]. Since the Baldwin-Barth turbulence model is used in the 
present work, the corrected angle of attack, a comp = 2.31°, suggested by Barth [7] 
is used. The far field and near field meshes are presented in Figures 7.8 and 7.9, 
respectively. The far field boundary is set at a radius of 5 chords from the airfoil. 
At this boundary the vortex far field boundary condition [83] is applied. A no-slip 
adiabatic wall condition and a zero normal pressure gradient are imposed at wall. 

The initial local structured mesh is constructed through the specification of a 
normal wall mesh scale (6„+ = Ay y+ = 8.6xl0 -6 chord) and a streamwise wall 
mesh scale ( N , — 25 and S w — A w <S y + = 0.0215 chord). In the far field the un- 
structured mesh is uniform with a length scale of = 0.5 chord. The refinement 
ratio equals 39 for the first adaptation and 2 for every additional refinement. The 
reason of using a large ratio in the first adaptation is to decrease the minimum 
length scale to the streamwise wall mesh scale so that the number of adaptation 
cycles can be reduced. Figures 7.9 to 7.13 present the resulting solution adapted 
meshes. Mesh statistics and force coefficient information for each mesh are listed 
in Table 7.3. Note this case takes about 2.5 hours of Cray-YMP CPU time. 


in r 
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Comparing the initial mesh solution with the final mesh solution as shown in 
Figures 7.14 and 7.15 shows improved shock wave resolution for the adapted solu- 
tion. Figures 7.16 to 7.21 show detailed mesh and flow structures in the shock wave 
and trailing edge regions. The present adaptive remeshing scheme does a good 
job in resolving the shock wave, the viscous wake and the boundary layer growth 
resulting from the strong shock wave. Figures 7.22 and 7.23 show the compar- 
isons of surface pressure coefficient with the experimental data for the initial and 
final meshes, respectively. These plots show a slight difference in the leading edge 
suction surface pressure distribution between the computed and the experimental 
results. This may be due to the fact that the transition trip was not simulated in 
the numerical solution. Further, the present solution does not accurately predict 
the shock wave location. This leads to the 3 % difference between the computed 
and measure lift coefficient (see Table 7.3). Comparing the initial mesh surface 
friction coefficient shown in Figure 7.24 with the final mesh surface friction coeffi- 
cient in Figure 7.25 shows a change of upper and lower surface friction coefficient 
as the mesh is refined. The change of upper surface friction coefficient is due to the 
resolution of the shock wave. A possible explanation to the change of lower surface 
friction coefficient is that the turbulence effect was not fully developed in the initial 
mesh solution yet. 

Since the flow conditions and turbulence model used here are the same as those 
used by Barth [7], the present solution shown in Figure 7.23 is compared to Barth’s 
solution shown in Figure 7.27. The present solution does not predict the same the 
shock wave location as Barth. In order to study the viscous effect on the shock 


Table 7.2 Flow conditions for the RAE2822 airfoil. 


A/oo 

^meai 

^comp 

(° K) 


LrrS (m) 

0.725 

2.92° 

2.31® 

293 

6.5x10® 

0.61 

Note: 

Experimental data is available in Reference [20]. 
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wave location, the inviscid solution for the same flow conditions, M«> = 0.725 and 
a = 2.31, is also computed. The inviscid surface pressure coefficient is presented 
as dashed line in Figure 7.26. The inviscid shock wave location is downstream of 
the experimental location as expected. Note that the artificial dissipation plays an 
important role in the accuracy of viscous flow solution as shown in the flat plate 
case. Since the difference between the shock wave locations of the expenment and 
viscous solution is 4 % of chord upstream, This shift may be due to too much level 
of artificial dissipation in the viscous flow region. However, comparing the present 
viscous results with the results of the 1987 Viscous Transonic Airfoil Workshop [32] 
as shown in Figures 7.23, 7.25, 7.28, and 7.29, the present results are well in the 
range of the workshop results. This indicates that this is very difficult test case to 

accurately predict. 

Figures 7.30 to 7.33 show comparisons of nondimensional velocity profiles at 
two stations on the airfoil upper surface. Since it is difficult to interpolate data 
within unstructured mesh regions, only the nondimensional velocity profiles within 
the local viscous mesh are shown. The nondimensional velocity is defined as q/U e , 
the ratio of the local velocity to the edge velocity. Because it is difficult to define 
an edge velocity on the unstructured mesh, the edge velocity is estimated using the 
procedure described in Reference [20], where the velocity ratio UJU^ is given as: 

U t _ M^ A -f 0.2M^ y /2 (7.2) 

Uoo Moo \ 1 + 0.2M* / 


with the edge Mach number Af e defined as: 


M e = 


( 1 + 0-2 Ml 

V VO+0.7 Ml(C p ) e ) 2/1 



(7.3) 


where Moo is the freestream Mach number and (C p ) e is the edge pressure coefficient. 
With the assumption of constant static pressure through the boundary layer, (C P ) e 
can be replaced with the surface pressure coefficient. The velocity ratio U e /Uoo is 
then evaluated using the above expressions. The plots show little difference between 
the initial and final mesh solutions. This is due to the fact that the normal wall 
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mesh scale is nearly constant during adaptation. Therefore, the resolution of viscous 
region is nearly constant during remeshing. 

Convergence histories for the pressure portion of the lift and drag coefficients 
are shown in Figures 7.34 and 7.35. Figures 7.36 and 7.37 present the convergence 
histories for the average residuals, \8{pu)/8t\ and \8Rtf8t\, respectively. The turbu- 
lence quantities take about 1,500 iterations to develop and the complete procedure 
converged in 8,000 iterations. These plots also show that all adapted mesh solutions 
do not satisfy the convergence criteria based on the average residual \8{pu)/8t\. 



Table 7.3 RAE2822 airfoil: = 0.725, a = 2.31°, Re L = 6.5x10® 

C lmt = 0.743, and C dm „ = 0.0127. 


MESH 


nodes 


elements 




Cb 


C 




8 — ./chord 


8 jomx / ^min 


Ay w «u(xl0 6 ) 


0 

1 

2 

3 

4 

6,568 

7,168 

8,401 

10,420 

14,348 

12,934 

14,133 

16,597 

20,626 

28,462 

0.7350 

0.7209 

0.7195 

0.7179 

0.7177 

0.0068 

0.0061 

0.0059 

0.0057 

0.0056 

0.0040 

0.0037 

0.0035 

0.0035 

0.0034 

0.00034 

0.00032 

0.00030 

0.00030 

0.00030 

0.5 

0.5 

0.5 

0.5 

0.5 

1 

39 

78 

156 

312 

4.28 

3.32 

3.69 

3.92 

4.02 

2 

2.5 

2.5 

2.5 

2.5 

2,500 

2,500 

2,500 

2,500 

2,500 

3,012 

1,932 

1,289 

874 

1,793 
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Figure 7.8 RAE2822 airfoil initial mesh: far field view. 
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Figure 7.9 RAE2822 airfoil initial mesh: 6,568 nodes and 12,934 elements. 
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Figure 7.10 RAE2822 airfoil 1 st adapted mesh: 7,168 nodes and 14,133 elements. 
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Figure 7.11 RAE2822 airfoil 2 nd adapted mesh: 8,401 nodes and 16,597 elements. 
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Figure 7.12 RAE2822 airfoil 3 rd adapted mesh: 10,420 nodes and 20,626 elements. 



Figure 7.13 RAE2822 airfoil 4 th (final) adapted mesh: 14,348 nodes and 28,462 

elements. 
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Figure 7.14 RAE2822 airfoil Mach number contours on the initial mesh: 
cmin = 0.0, cmax = 1.20, and inc = 0.05. 


Figure 7.15 RAE2822 airfoil Mach number contours on the final mesh 
cmin = 0.0, cmax = 1.20, and inc = 0.05. 
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Figure 7.16 RAE2822 airfoil: shock region blowup of the final mesh. 



Figure 7.17 RAE2822 airfoil: trailing edge blowup of the final mesh. 




Figure 7 18 RAE2822 airfoil Mach number contours (cmin = 0.0, cmax 
and inc = 0.05): final mesh shock region. 


cure 7.19 RAE2822 airfoil Mach number contours (cmin - 0.0, cmax 
and inc = 0.05): final mesh trailing edge region. 
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Figure 7.20 RAE2822 airfoil velocity vector: final mesh shock region. 



Figure 7.21 RAE2822 airfoil velocity vector: final mesh trailing edge region 




Figure 7.23 RAE2822 airfoil surface pressure coefficient for the final mesh: solid 
line - numerical solution and symbol - experimental data. 
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Figure 7.24 RAE2822 airfoil surface friction coefficient for the initial mesh: solid 
line - numerical solution at lower surface, dashed line - numerical solution at 
upper surface and symbol - experimental data at upper surface. 



*/c 


Figure 7.25 RAE2822 airfoil surface friction coefficient for the final mesh: solid 
line - numerical solution at lower surface, dashed line - numerical solution at 
upper surface and symbol - experimental data at upper surface. 
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Figure 7.26 RAE2822 airfoil surface pressure coefficient for the final mesh: solid 
line - numerical turbulent solution, dashed line - numerical inviscid solution, and 

symbol - experimental data. 


U 0.0- 



Figure 7.27 RAE2822 airfoil surface pressure coefficient: Barth’s numerical 
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Figure 7.28 RAE2822 airfoil surface pressure coefficient: the 1987 viscous 
transonic airfoil workshop compendium of results [32]. 



Figure 7.29 RAE2822 airfoil surface friction coefficient: the 1987 viscous 
transonic airfoil workshop compendium of results [32]. 
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Figure 7.30 RAE2822 airfoil nondimensional velocity profile at x=0.319 for the 
initial mesh: solid line - computed results and symbol - experimental data. 



Figure 7.31 RAE2822 airfoil nondimensional velocity profile at x=0.319 for the 
final mesh: solid line - computed results and symbol - experimental data. 
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7.3 NLR Two- Element Airfoil 

The NLR two-element airfoil case illustrates the ease with which viscous flow 
problems on arbitrary multiply-connected regions can be solved. This two-element 
airfoil is composed of a NLR7301 main airfoil with a trailing edge flap. Experiments 
have been done for various angles of attack with a single flap deflection angle of 20° 
by van den Berg [84]. For the case presented the flow conditions are summarized 
in Table 7.4. This flow has a rapid expansion around the leading edge. There is 
a small laminar separation bubble on the upper surface starting from 2.6 to 4 % 
chord according to flow visualizations [84]. For such a small separation bubble, it 
is very important to have correct level of turbulent viscosity to resolve it. 

The far field and near field meshes are shown in Figures 7.38 and 7.39, respec- 
tively. The far field boundary is placed at a radius of 5 chords from the airfoil. 
The vortex far field boundary condition [83] is specified at this boundary. Along 
wall surfaces, no-slip condition with adiabatic wall and zero normal pressure gra- 
dient boundary conditions are applied. The initial C-type local structured mesh is 
generated around the airfoil through the specification of a normal wall mesh scale, 
8 y + = 2xl0 -s chord, with a streamwise wall mesh scale (N, = 40 and 8 W = 0.015 
chord). In the far field the unstructured mesh is uniform with a mesh scale of 
= 0.5 chord. The sequence of adapted meshes are shown in Figures 7.39 to 
7.41. Mesh properties and force coefficients are listed in Table 7.5. Comparing 
the final mesh Mach number contours with the initial mesh solution as shown in 
Figures 7.42 and 7.43 shows improved resolution of the boundary layer edges. The 
mesh and flow at the leading edge and in the airfoil-flap gap region are shown in 


Table 7.4 Flow conditions for the NLR two-element airfoil. 


Moo 

a 

BBS 

Flap gap % c 

Too C K) 

Rei 

Lre f (m) 

0.185 

6.0° 

20.0° 

2.6 

293 

2.51x10® 

0.57 

Note: 

Experimental data is available in Reference [84]. 



[If I 

















169 


Figures 7.44 to 7.47. Although these plots show an accurate prediction of stag- 
nation point location, the leading edge separation bubble is not resolved. This is 
also confirmed by the positive value of the friction coefficient distribution along the 
surface as shown in Figure 7.53. Since the refinement parameter based on surface 
curvature does a good job in resolving the leading edge region, plots of the initial 
and final mesh surface pressure distribution are almost exactly the same (see Fig- 
ures 7.50 and 7.51). The surface friction coefficient distribution does improved as 
the mesh is refined as shown in Figures 7.52 and 7.53. Similar improvements are 
shown in the velocity plots (see Figures 7.54 to 7.57). This improvement is due to 
a better resolution of near-wall viscous regions since a smaller wall distance Ay„«n 
is obtained as the solution develops (see Table 7.5). 

Figures 7.58 and 7.59 present the convergence histories for the pressure portion 
of the lift and drag coefficients. Convergence histories for the average residual 
of \8{pu)/8t\ and \8Rt/8t\ are given in Figures 7.60 and 7.61, respectively. This 
case takes about 5.7 hours of Cray-YMP CPU time. Although the difference in 
the force coefficient between last two meshes does not fall within the convergence 
criteria (.002), the adaptation cycle has been stopped due to the huge computation 
time required for continued adaptations. Further refinement of the mesh will not 
improve the solution since the wall distance has converged to a constant value on 
the last two adaptations (see Table 7.5). The computed total lifting coefficient is 
about 1.25 % larger than the measured value. 
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Table 7.5 NLR two-element airfoil: M = 0.185, a = 6.0°, 6fl*p = 20°, 
Re L = 2.51X10 6 , C tmw = 2.416, and C dtn ^ = 0.0229. 


MESH 

0 

1 

2 

nodes 

9,510 

11,579 

13,666 

elements 

18,708 

22,793 

26,931 


2.3983 

2.4345 

2.4415 

c dp 

0.0470 

0.0411 

0.0341 

Cl, 

0.00358 

0.00415 

0.00491 

C*. 

0.00073 

0.00084 

0.00094 

/chord 

0.5 

0.5 

0.5 

^mAx/^min 

1 

39 

78 

Ay w «u(xl0 6 ) 

10.0 

5.56 

4.81 

2/pr 

2 

2.5 

2.5 


750 

750 

750 

Cray-YMP CPU (sec) 

9,770 

3,362 

7,411 


Ill 1 























171 



Figure 7.38 


NLR two-element airfoil initial mesh: far field view. 
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Figure 7.39 NLR two-element airfoil initial mesh: 9,510 nodes and 18,708 

elements. 
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Figure 7.40 NLR two-element airfoil l* 4 adapted mesh: 11,579 nodes and 22,793 

elements. 
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Figure 7.41 NLR two-element airfoil 2 nd (final) adapted mesh: 13,666 nodes and 

26,931 elements. 
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Figure 7.44 NLR two-element airfoil: leading edge blowup of the final mesh. 



Figure 7.45 NLR two-element airfoil: airfoil-flap gap blowup of the final mesh. 
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Figure 7.50 NLR two-element airfoil surface pressure coefficient for the initial 
mesh: solid line - numerical solution and symbol - experimental data. 



Figure 7.51 NLR two-element airfoil surface pressure coefficient for the final 
mesh: solid line - numerical solution and symbol - experimental data. 
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-0.0080-1 

- 0.20 


IN 
— \ 


1 

i 



Figure 7.52 NLR two-element airfoil surface friction coefficient for the initial 
mesh: solid line - numerical solution at lower surface, dashed line - numerical 
solution at upper surface and symbol - experimental data at upper surface. 
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Cf 
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- 0.20 


Figure 7.53 NLR two-element airfoil surface friction coefficient for the final mesh: 
solid line - numerical solution at lower surface, dashed line - numerical solution at 
upper surface and symbol experimental data at upper surface. 
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Figure 7 54 NLR two-element airfoil nondimensional velocity profile at x-0.88 for 
the initial mesh: solid line - computed results and symbol - experimental data. 



Figure 7.55 NLR two-element airfoil nondimensional velocity profile at x-0.88 for 
the final mesh: solid line - computed results and symbol - experimental data. 
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Figure 7.56 NLR two-element airfoil nondimensional velocity profile at x= 1.11 for 
the initial mesh: solid line - computed results and symbol - experimental data. 



Figure 7.57 NLR two-element airfoil nondimensional velocity profile at x=l.ll for 
the final mesh: solid line - computed results and symbol - experimental data. 









7.4 VKI LS82-1 Turbine Cascade 


This turbine cascade case is 


referred to as case 1 in the 1982 VKI Work- 


shop [77], Experimental results for the surface isentropic Mach number distribution 
and Schlieren photographs for different exit pressures are given by Sieverding, and 
Van Hove et al. [77]. The isentropic Mach number is defined as 


M» = 


2 

1-1 



(7.4) 


The blade characteristics and the flow conditions for the results presented are sum- 
marized in Table 7.6. The exit pressure is measured in an area averaged fashion 
and is given in terms of the isentropic exit Mach number. The ratio of the area av- 
eraged exit static pressure to the inlet total pressure is 0.3012 for the present value 
of (Af„w. This flow condition is of particular interest because there is a strong 
oblique shock wave propagating to downstream boundary as shown in Schlieren 
photograph (Figure 7.68a). The total exit velocity is supersonic with a subsonic 
axial component. Therefore, one boundary condition is stiff needed at the exit 
boundary. In addition, the flow in this case is choked, making the boundary layer 
resolution very important. Note that experimental results are measured using a 
finite number of blades with tailboard placed on the upper boundary as shown in 
the Schlieren photograph. The use of tailboards results in reflected shock waves 
which propagate to the blade row. These reflected shock waves are not presented 

in an infinite cascade problem. 


Table 7.6 Flow conditions for the VKI LS82-1 turbine cascade. 


Min 
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An 
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Pitch/ chord 

0.085 

1.43 



278.0 

2.51x10* 

0.78 


Note: • Experimental data is available in Reference [77]. 


• Nonuniform pressure present at the exit. 
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The initial mesh is given in Figure 7.62. The upstream and downstream mesh 
boundaries are placed one half an axial chord from the blade leading and trailing 
edges, respectively. At the upstream boundary total pressure, total temperature, 
and absolute flow angle are specified as boundary conditions. At the downstream 
boundary the area averaged exit static pressure is imposed as described in Chapter 
3. No-slip condition, adiabatic wall, and zero normal pressure gradient are enforced 
on the blade and periodicity is imposed at the upper and lower boundaries. Spec- 
ifying a normal wall mesh scale (£„+ = 6xl0 -5 chord) and a streamwise wall mesh 
scale (N, = 8 and 6 W = 0.018 chord) produces the initial O-type local structured 
mesh around the blade. In the far field the unstructured mesh is uniform with a 
length scale of = 0.0385 (pitch). 

The sequence of solution adapted meshes are shown in Figures 7.62 to 7.65. 
The characteristics of each mesh axe tabulated in Table 7.7. Figures 7.66 and 7.67 
present the initial and final mesh solutions, respectively. The latter plot shows a 
good resolution of the oblique shock wave and viscous wake. Since the Schlieren 
photograph is based on the reflection of light passing through the flow, the density 
contours of the computed results are compared to the Schlieren photograph. Com- 
paring the final mesh density contours shown in Figure 7.68b with the Schlieren 
photograph in Figure 7.68a shows an accurate prediction of the strong oblique shock 
wave location and orientation, as well as the weak passage shock wave. Detailed 
mesh and Mach contours at the blade trailing edge are presented in Figures 7.69 
and 7.70. These plots demonstrate how the solution adaptive scheme resolves com- 
plex flow structures. Comparing the initial mesh surface Mach number distribution 
shown in Figure 7.71 with the final mesh solution in Figure 7.72 shows an improved 
solution as the mesh is refined. Figure 7.72 shows a slight difference between the 
numerical solution and experimental data on the uncovered portion of the blade 
suction surface. This is due to the interference of the reflected shock waves from 
the tailboard with the suction surface of blade 6 in the experimental data (see 
Figure 7.68a). 
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Although the solution converges after three mesh adaptations, further refine- 
ment is forced to improve the resolution of the weak shock wave. The 4 th adapted 
mesh is fine enough to resolve the vortex shedding at the blade trailing edge as 
shown in Figure 7.77. Since vortices are shed from the blunt trailing edge, this flow 
is unsteady. In order to obtain a realistic unsteady solution, the solution should be 
marched time accurately. This is too costly to be practical. In the present work 
a "pseudo” time averaged quantity is used to estimate the difference in loading 
coefficients between steady state and unsteady solutions. The term "pseudo” is 
used to distinguish it from the actual time average. Due to the use of local time 
stepping and implicit residual smoothing, the present calculation does not produce 
a time-accurate solution. The mesh in the unsteady wake flow region as shown in 
Figure 7.75 is quite uniform, therefore, the solution in this region may be nearly 
time-accurate. A "pseudo” time averaged loading coefficient can be estimated from 
the convergence histories of the solution on the 4** adapted mesh. It is believed 
that the trailing edge shedding has a small effect on the loading coefficient. Since 
the loading coefficient is a major concern in practical applications it might be rea- 
sonable to accept the steady state solution. This is confirmed by the convergence 
histories for the axial and normal force coefficients shown in Figures 7.81 and 7.82, 
respectively. The fluctuating part of the loading coefficient is less than 1 % of the 
steady state or "pseudo” time averaged value. Comparing the 4 th adapted mesh 
solutions as shown in Figures 7.79 and 7.78 with the steady state solution in Fig- 
ures 7.72 and 7.78 also shows little difference in the two solutions. The total CPU 
time for three mesh adaptations is about 2.2 hours of Cray-YMP CPU time. 
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Figure 7.62 VKI LS82-1 turbine cascade initial mesh: 4,610 nodes and 8,780 

elements. 
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Figure 7.63 VKI LS82-1 turbine cascade 1** adapted mesh: 5,212 nodes and 9,918 

elements. 
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Figure 7.64 VKI LS82-1 turbine cascade 2” 1 adapted mesh: 7,783 nodes and 

14,943 elements. 
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Figure 7.66 VKI LS82-1 turbine cascade Mach number contours on the initial 
mesh: cmin = 0.0, cmax = 1.8, and inc = 0.05. 
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Figure 7.67 VKI LS82-1 turbine cascade Mach number contours on the 3 rd 
adapted mesh: cmin = 0.0, cmax = 1.8, and inc = 0.05. 
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(a) (b) 


Figure 7 68 VKI LS82-1 turbine cascade: (a) Schlieren photographs (Sieverdmg 
[77]). (b) Density contours on the 3 rd adapted mesh, cmin = 0.17, cmax - 0.99, 

and inc = 0.02. 
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Figure 7.69 VKI LS82-1 turbine cascade: trailing edge blowup of the 3 rd adapted 

mesh. 
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Figure 7.71 VKI LS82-1 turbine cascade surface isentropic Mach number for the 
initial mesh: solid line - numerical solution and symbol - experimental data. 



Figure 7.72 VKI LS82-1 turbine cascade surface isentropic Mach number for the 
3 rd adapted mesh: solid line - numerical solution and symbol - experimental data. 
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■Figure 7.77 VKI LS82-1 turbine cascade velocity vector: 4 th adapted mesh 

trawling edge region. 
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Figure 7.78 VKI LS82-1 turbine cascade Mach number contours on the 4 th 
adapted mesh: cmin = 0.0, cmax = 1.80, and inc = 0.05. 


Ill I 



205 



Figure 7.79 VKI LS82-1 turbine cascade surface isentropic Mach number for the 
4 th adapted mesh: solid line - numerical solution and symbol - experimental data. 



Figure 7.80 VKI LS82-1 turbine cascade surface friction coefficient for the 4 th 

adapted mesh. 
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7.5 Allison Tandem Blade Cascade 

This test case demonstrates the present solution adaptive scheme for a complex 
blade design. There is no experimental data available. Since the flow is three- 
dimensional, through-flow information is required as input to the present quasi- 
three-dimensional analysis. The streamsurface location and thickness data used 
in the present calculations were obtained from Reference [60]. The adapted mesh 
solutions presented here are for 70 % span station measured from the shroud. This is 
the same case for which the inviscid solution of section 6.6 was discussed. The inflow 
and outflow conditions are determined from the through-flow analysis. The exit 
static pressure is adjusted until a correct incident is obtained. The flow conditions 
for the 70 % span station are summarized in Table 6.12. 

The initial mesh is given in Figure 7.85. The upstream and downstream mesh 
boundaries are located one axial chord from the blade edges. At the upstream 
boundary total pressure, total temperature, and absolute flow angle are imposed. 
At the downstream boundary the exit static pressure is specified. At blade surfaces 
a no-slip condition with adiabatic wall and zero normal pressure gradient is imposed. 
The initial O-type local structured mesh is constructed using a normal wall mesh 
scale 6 V+ = 9.48X10" 5 chord and a streamwise wall mesh scale (N. = 20 and 
S w = 0.055 chord). In the far field the unstructured mesh is uniform with a length 

scale of £m&x ~ 0.046 chord. 

The initial mesh solution is massively separated on the suction surface start- 
ing from 80 % of chord of the first blade to all the way downstream as shown in 
Figure 7.87. Similar separated flow solutions have been obtained by another inves- 
tigator [24]. Comparing the turbulent flow results shown in Figure 7.87 with the 


Table 7.8 Allison tandem blade cascade far field flow conditions. 
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inviscid flow results in Figure 6.61 shows the importance of modeling the viscous 
effects. 

It is meaningless to perform a further mesh adaptation on this kind of flow 
structure. The mesh adaptation is performed one time to illustrate how the solution 
adaptive scheme resolves such separated flow structures. Figure 7.94 shows mesh 
refinment on the separated and wake flow regions. The present mesh refinement 
parameters do a good job in resolving viscous flow structures. The mesh statistics 
and force coefficient information is summarized in Table 7.9. This case demonstrates 
that the unstructured flow solver can easily produce solutions for flow problems 
which are difficult to predict using structured mesh approach. 



210 


Table 7.9 Allison tandem blade cascade 70 % span station: A/ ex it = 0.71, 

= 0.0°, and Re L = 5.3x10 s . 
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Figure 7.85 Allison tandem blade cascade 70 % span from shroud initial mesh: 

far field view. 



Figure 7.86 Allison tandem blade cascade 70 % span from shroud initial mesh: 

5,417 nodes and 10,481 elements. 
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Figure 7.87 Allison tandem blade cascade 70 % span from shroud Mach number 
contours on the initial mesh: cmin =0.0, cmax= 1.45, and inc = 0.05. 


Figure 7.88 Allison tandem blade cascade 70 % span from shroud: leading edge 

blowup of the initial mesh. 



Figure 7.89 Allison tandem blade cascade 70 % span from shroud Mach number 
contours (cmin = 0.0, cmax = 1.45, inc = 0.05): initial mesh leading edge region. 


Figure 7.90 Allison tandem blade cascade 70 % span from shroud: tandem blade 

gap blowup of the initial mesh. 



Figure 7.91 Allison tandem blade cascade 70 % span from shroud Mach number 
contours (cmin = 0.0, cmax = 1.45, inc = 0.05): initial mesh tandem blade gap 

region. 
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Figure 7.92 Allison tandem blade cascade 70 % span from shroud velocity vector: 

initial mesh leading edge region. 



Figure 7.93 Allison tandem blade cascade 70 % span from shroud velocity vector: 

initial mesh tandem blade gap region. 
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Figure 7.94 Allison tandem blade cascade 70 % span from shroud 1** adapted 
mesh: 9,367 nodes and 18,186 elements. 
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8. SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS 

A general solution adaptive unstructured remesh scheme has been implemented 
and verified for a wide range of both inviscid and viscous two-dimensional and 
quasi-three-dimensional turbomachinery flows. The solution adaptive scheme used 
here incorporates an explicit finite- volume Runge-Kutta time-marching flow solver 
and an advancing front mesh generation scheme. The mesh is adapted by periodic 
remeshing as the solution evolves. In the present approach, an edge-based local 
coordinate system has been introduced for the formulation of the stability criteria 
and eigenvalue scaling of the artificial dissipation. A general stability analysis of an 
unstructured scheme is not possible due to the unstructured mesh connectivity. By 
employing the edge-based coordinate system, an analogy to the structured formu- 
lations is obtained. To enhance the convergence rate to steady state, an improved 
iterative method for solving the implicit residual averaging operator is formulated 
which allows larger local time steps to be used than those obtained from standard 
Jacobi iterative method. 

In the solution adaptive scheme, new mesh refinement parameters have been 
derived based on a combination of the local surface curvature and the gradient of 
flow speed. The coupling of geometric and flow field information results in an ac- 
curate and efficient adaptation criterion for problems with complex flow structures 
and complex geometries. The present work uses a local structured triangular mesh 
within viscous flow regions. The solution adaptive remesh scheme is general and is 
not restricted to this local structured mesh. 


hit: 
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The major contributions of the present research are summarized below. 

• A solution adaptive unstructured remesh algorithm has been developed for 
turbomachinery applications. 

• New mesh refinement parameters have been formulated for an accurate and 
efficient resolution of complex flow problems. 

• Improvements to the basic cell- vertex finite- volume Runge-Kutta scheme have 
been made in the following areas: 

- Introduction of an edge-based local coordinate system for the formula- 
tion of the stability criteria and eigenvalue scaling artificial dissipation. 

- Development of an improved iterative method for the implicit residual 
averaging equation. 

• Demonstration of the excellent performance of the solution adaptive unstruc- 
tured remesh scheme. 

8.1 Summary of Results 

The solution adaptive remesh scheme has been demonstrated with the solu- 
tion of various two-dimensional and quasi-three-dimensional inviscid and turbulent 
flows. Each of these cases has different flow features which need to be resolved for 
an accurate solution. In the cases presented, the adaptive remesh scheme was used 
to automatically resolve the flow features in the region of shock waves, expansion 
waves, stagnation regions, viscous layers and viscous wakes. These cases illustrates 
the automatic detection of the location, the size and the orientation of the impor- 
tant flow features through the use of the solution adaptive remesh algorithm. The 
minimum mesh scale of the sequence of solution adapted meshes is user-specified 
and can be defined in an efficient fashion to enhance the convergence of mesh adap- 
tation cycles. This leads to a great advantage over the mesh refinement approach 
where mesh cells are subdivided in a constant ratio of 2. 
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Flow problems of arbitrary multiply-connected regions are no more difficult 
to solve than flow over an isolated body as shown by multi-element airfoil and 
tandem blade cascade cases. Small finite radius leading and trailing edge regions 
are also resolved as shown in the NASA Rotor 67 30%-span station case (Figure 
6.51). These cases clearly illustrate the great flexibility of the solution adaptive 
unstructured remesh scheme for geometrically complex flow problems. 

Even when complex flow structures develop within the solution domain, the 
solution adaptive remesh algorithm accurately resolves them as shown in the Den- 
ton supersonic wedge cascade and the VK1 turbine cascade cases. These complex 
flow structures are very difficult to be accurately resolved using structured mesh 
approaches. Comparing the final solutions with the corresponding analytical or 
experimental data has shown that the important flow structures are accurately 
resolved without any prior knowledge. There is also a great saving in the computa- 
tional effort and storage requirements since mesh points are not wasted in regions 
where they are not required. Furthermore, any level of the accuracy of solutions 

can be achieved as the mesh continues to be refined. 

The effect of artificial dissipation on the solution of turbulent flows has also 
been studied. While the adaptive artificial dissipation works well in the solution 
of inviscid flows, it causes unacceptable errors in the solution of turbulent flows as 
shown in the flat plate case. Shutting off artificial dissipation within boundary layer 
regions produces an accurate prediction of the turbulent flow, but it is impossible 
to turn off the artificial dissipation completely for general flow problems. With 
proper scaling of the artificial dissipation, it has been shown that these errors can 

be greatly reduced. 

In a preliminary study of directional mesh adaptations, it has been found that 
the advancing front method produces highly distorted meshes when mesh cell as- 
pect ratios are greater than 20. The aspect ratio limitation in the advancing front 
method limits its use in viscous flow regions since high aspect ratio meshes are re- 
quired for efficient viscous flow calculations. In the present work, a local structured 
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mesh is used in viscous regions. Adopting a local structured mesh in viscous flow 
regions yields several advantages over a simple unstructured mesh approach. First, 
high aspect ratio mesh cells are easily achieved. This in turn yields an efficient 
mesh point dis tribution for the solution of turbulent flow problems. Second, vis- 
cous mesh scale is easily controlled, which reduces the possible stability problem 
caused by insufficient resolution of the near-wall region for turbulent flows. 

For turbomachinery cascade calculations, repeated mesh refinement resolves 
vortex shedding at blunt trailing edges as shown in the VKI turbine cascade case. 
Since vortex shedding is an unsteady phenomena, a time-accurate solution scheme 
is required to obtain the correct solution. However, the present scheme can be 
marched time- accurately, this requires an enormous computational effort. Since 
the fluctuation in force coefficients is less than 1% of the steady state or pseudo 
time averaged force coefficients, the steady state solutions computed for this case 
are considered a good engineering approximation to the real flow. 

The present solution adaptive scheme has been shown to accurately predict 
complex turbulent flows, although the explicit time-marching scheme requires a 
great amount of computational time for solving turbulent flow problems due to 
the requirement of fine meshes in the near-wall regions. Nevertheless, the results 
demonstrate the capability of the present scheme in analyzing the complex turboma- 
chinery problems which can not be solved using classic structured mesh approaches 
(e.g., plate splitter, tandem blade, etc.). 

8.2 Future Work 

This thesis has successfully demonstrated a flexible solution adaptive remesh 
scheme for the accurate solution of complex flows. Further research in the following 
areas is suggested to enhance the efficiency of the present approach: 

-1. Multi-grid acceleration: 

Multi-grid acceleration should be incorporated into the present solver. Multi- 
grid methods have been demonstrated for unstructured triangular mesh Euler 
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and Navier-Stokes solutions by Mavriplis [54]. While multi-grid has been 
found to be more computationally intensive for triangular mesh formulations 
due to the the complexity of interpolation and projection operations, there is 
still a significant reduction in the overall computational work. 

2. Alternate flow solvers: 

For viscous flow calculations the computation time becomes extremely large 
due to the stability criteria of the present Runge-Kutta scheme. This may be 
overcome through the adoption of an implicit formulation. There has been 
some promising work in this area by Batina [10], Barth [7], and others, but 
much work remains to be done. 

3. Directional adaptation: 

Since many flow structures, such as viscous layers and shock waves, involve 
multiple length scales where the change in flow properties along one direc- 
tion is much larger than the change in other directions, it is inefficient to 
use equilateral triangles over the domain. A directional adaptation method, 
which allows high aspect ratio triangles to be aligned with the flow features, 
will further reduce the computation expense and storage requirements en- 
countered in these regions. The present mesh generation scheme is capable of 
generating directionally adapted mesh with cell aspect ratios as high as 20. 
With further development of the refinement parameters this could be incor- 
porated in inviscid flow regions. The high aspect ratio meshes required for 
efficient viscous flow calculations require a new unstructured mesh generation 
scheme. 
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Appendix A: Derivation of the Quasi-Three-Dimensional Navier-Stokes Equations 


The vector form of the three-dimensional Navier-Stokes equations may be ex- 
pressed as follows. 


! + v,?=° 


El + sz = l v .i 

Dt P P 


and 


dpE 


(A.1) 


(A.2) 


(A.3). 


_ + V • pVH = V • («. V T + i ■ V) 
at 

The above equations represent the physical conservation law of mass, momentum 
and energy for a fluid flow without body force and external heat transfer. For 
analysis of general flow problem it is useful to consider flow on the generalized 
curvilinear coordinate system because they can be transformed to any coordinate 

system. 

Let (xi,x 2 ,x 3 ) be a set of generalized curvilinear coordinates and (ei,e 2 ,e 3 ) 
be the set of corresponding unit vectors. The expressions for the vector operators 
appearing in the Navier-Stokes equations are described below. 

The gradient of a scalar <f> is described by 


, 1 d<J> , 1 d<f> , 1 d<j>^ 

V <f> = — e l + — "5— e 2 + — e 3 


k\ dx\ A ' ^2 9x2 

The gradient of a vector A is given by 


/i 3 dxz 


V A = [An] 

where Aij is the ij th element of a tensor matrix A and is described as 

1 dA± c A k 9hj Aj dh : 

** hj dxj IJ kz=zl h{hk dxk hihj dx{ 

where Sij is the Kronecker delta function. 
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Applying Eq.(A.S) 

(A-V)B = 


ej t dB\ dB 2 dB 3 ^ A 2f dh 2 B 2 dhiBi x 

T (^ 1 - 3 — + A 2 -^ — + A 3 - — ) - — (— ) 

h\ OX\ UX\ QX\ U2 UX\ UX 2 
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The divergence of a vector A is given by 

V *A = 7— J— 7- p — (^2^3 Ai) + — — (A3A1A2) + — (^1/12^3) 

A1/I2A3 [< 7 Xi CrX2 C/X3 

The divergence of the stress tensor a is expressed as 

V * = e, {vhl + A-(M 3 <t, 2 ) + A-fA.A,. 
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The stress tensor can be described as 

ff = /f(^ + f)+Mv -v 

where S denotes the velocity gradient tensor matrix. 


(A. 10 ) 
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The metric coefficients h\, & 2 , and h$ are defined as 



(A.ll) 

In order to obtain the Navier-Stokes equations for the quasi-3D blade to blade 
flow, we assume the coordinate system (m,0) is on the streamsurface and the third 
component n is perpendicular to the plane (m, 0). With this third coordinate, we 
can define a three-dimensional curvilinear coordinate system using (m, 6, n). In this 
coordinate system we have V • e n = 0 or V n = 0, so the velocity may be described 
as 

v = (v m ,v a ,o) (A.i2) 

If we let the streamtube height h be defined in the direction of n and measured in 
the units of the variable n, then the metric coefficients may be obtained as 

A! = l ^ l=r ’ * 3 = l fi; l = '‘ (A - 13) 

Substituting the vector operators, metric coefficients, and velocity into the vector 
form of the Navier-Stokes equations gives the quasi-3D Navier-Stokes equations as 
follows. 

Continuity equation: 
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(A.14) 
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m- momentum equation: 
dV, 
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Equations (A.15) and (A.16) can be cast in the conservation form using the conti- 
nuity equation(A.14). 

Performing a linear combination of Eq.(A.15)-rA/> + Eq.(A.14)*Kn, the m- 
momentum equation is obtained as 

J^(r/i pV m ) + ^[rfe(pV£ + p)] + Qg(hpV 9 V m ) 

' x 1 dr / . 1 dh 

(pV} +p- + (p - " 33 ) 5 -^ 

(A.18) 

Similarly, a linear combination of r- (Eq.(A.16)-rhp + Eq.(A.14)-V,) results in 
the conservation form of the 6 - momentum equation. 

r I + + gff^ pV ^ + + ^ p ^ mVe dm } 

= T ~L. {rk ° n) + (r '“ T,j) ^ + r m {h ' rr,) 

-> §j(rhpV,r) + A( rW V„r) + ^ [A(pV. a + p)r] 

= ^(r/t< 7 12 r) + (A. 19) 

Therefore, the quasi-3D Navier-Stokes equation may be expressed as 

djr hU) d(rhF) d(hG) _ d{r hR) 9(hg) rh g (A 20) 

dt + dm 86 dm 80 

where 
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Employing Eq.(A.lO) the stress terms axe 

o n = 2nd m V m + \\ 7 -V 

o 22 = 2ft(d e V e + V m -^)/r + XvV 

033 = 2 ^( 1 ^) + A V-V 

an = fi(d m V 0 -V e (~) + -deV m ) (A.22) 

r am r 

We may also write the equations in the relative frame which the coordinate system 
rotates with blade rows. Let fi be the rotating speed of blade rows, 6 be the relative 
angle with respect to the blade rows, and W$ denote the relative tangential velocity, 
then we have 

6 = 6-nt, W 8 = V e - rft 


and 


d__d_ dd&_ = d__ Q d_ 
dt~ dt + d0dt dt U d0 


Substituting the above relation into equation (A.20), we have 
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Appendix B: Quasi-Three-Dimensional Favre- Averaged Navier-Stokes Equations 


In compressible flows, the density- weighted average suggested by Favre [27 , 28] is 
helpful to simplify the formulation of the mean flow equations. Favre averaging is a 
hybrid average method which uses density-weighted averaging on all fluid properties 
except for pressure and density on which an ensemble (or time) averaging is used. 
An ensemble-averaged quantity, represented by angular brackets <>, is defined as 

</>=4e/. < b » 

1=1 


and a density- weighted averaged quantity, represented by a tilde, is defined as 

(Pf) 


f = 


(P) 


(B.2) 


where f is some indenpendent variable. Follow the Reynolds averaging procedure, 
the quantity is decomposed into the averaged and fluctuating parts. 


/ = {/)+/', and(/)=0 (B.3) 


and 

/ = / + /", and (pf) =0 but r ± 0 (B.4) 

Employing Favre averaging the vector form of the mean flow equations in the gen- 
eralized curvilinear coordinate system are [74] 

^ + V-((/>)^)=0 (B-5) 


<B6) 

and 

^ + V ■ = V • (*. V r + + ?•(»- {pV-V-i ) ) (B.7) 

-The mean flow equations are in a similar form of the full Navier-Stokes equations 
Eqs.(A.l)- (A. 3) except for the appearance of the unknown terms in the mean mo- 
mentum and energy equations. Applying the same procedures given in Appendix A, 
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the Favre-averaged N-S equations in the generalized curvilinear coordinate system 
can be converted to the desired quasi-3D relative coordinate system. 

In the particular (m, 0 , n) coordinate system, the mean-flow velocity V and the 
fluctuating velocity V" are given by 

V = V m e m + Vge 9 

V" = V - 9 = V m "e m + V e "e e + V n "e n (B.8) 


In the above expressions, 0 now represents the relative angle with respect to the 
blade rows. Using the coordinate transformation procedure given in the Appendix 
A, the mean-flow equations Eqs.(B.5)-(B.7) are converted to 
Mean-flow continuity equation: 

|W/>)) + ^( r M^) + ^(AW^)=0 (B-9) 

Mean flow m- momentum equation: 
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Mean-flow momentum equation: 
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These mean-flow equations may be rewritten in a condensed vector form as 
d(rh{U)) d(rh(F)) d(h(G)) (d( rh(R)) d(h(S)) \ __ LI ^ 13) 

aT" + am + ae \ a™ * ao ) ' ' ' . 
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